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ABSTRACT 


This investigation deals with the use of a irovable 
mass control system to stabilize an arbitrarily tumbling 
asymmetric vehicle about the maximum inertia axis. A 
first-order gradient optimization technique is used to 
minimize angular velocity components along the intermediate 
and minimum inertia axes . This method permits a wide 
range of initial guesses for mass position history. Mo- 
tion of the control mass is along a linear track fixed in 
the vehicle. The control variable is taken as mass 
acceleration with respect to body coordinates. Motion is 
limited to defined quantities and a penalty function is 
used to insure a given range of positions. Numerical 
solutions of the optimization equations verify that 
minimum time detumbling is achieved with the largest 
permissible movable mass, length of linear track, and 
positions of the mass on the two coordinates perpendicular 
to the linear motion. Also, the mass should oscillate, 
about the zero point, on an axis parallel to the major 
principal axis. A minimum mass solution is obtained 
by fixing the time at the largest feasible value. The 
optimal method permits detumbling in about one-fourth 
the time when compared to a force control law formulation 
available in the literature. Since stabilization may 
require hours, this reduction in time is very significant. 


In regard to minimum mass, the optimization permits the 
use of a much smaller mass for detumbling in the same 
time. This mass reduction is quite substantial since 
very large masses are required. Use of this control 
system for actual operations in space is feasible since 
the velocity and acceleration of the mass , and the power 
requirement, are low. It should be noted that the control 
technique utilizes an open loop solution in real time. 

In addition, the technique need not be restricted to 
attaining simple spin about the maximum inertia axis ; 
geometric axes may be specified. 
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NOMENCLATURE 


a = Inertial acceleration of the origin of 
coordinates which is fixed at the center 
of mass of the main body 

i: = Force on the contx*ol mass 

m 

= x component of J 
f = n column vector composed of f^ 

f ^ = Functions to which x s ^ are equal 
f - Resultant of external forces 

Fp = Function on which an inequality constraint 
is placed 

Fp max = Maximum permitted value of 

Fpfi xe d = Actual Fp max used in computation 

5 - Angular momentum of the system with respect 
to the origin of coordinates . 

§ = Angular momentum of the movable control 

a mass with respect to the origin of 
coordinates 

Hg = Angular momentum of the main body with 
respect to the origin of coordinates 

HB x »Hg 9 Hjj = Component of Sp along the x 5 y 9 and z 
y axes 9 respectively 

8 = Angular momentum of the system with 

respect to its own center of mass 

^ 

i s 3 ,k = Orthogonal unit vectors of coordinate frams 
fixed in the main body 

I ,1 ,1 = Moments of inertia of the main body 

x y z 

I,,. = Products of inertia of the main body 

xy xz y z 

I = Inertia dyadic of the main body 

I = Maximum moment of inertia of the system 

max J 
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NOMENCLATURE (Continued) 


I . = Minimum moment of inertia of the system 

mm 

I . . = (q x q) matrix required in the first-order 
^ gradient method 

I |T = q column vector required in the first-order 
^ gradient method 

I TT = Scalar required in the first-order gradient 
JJ method 

J - Performance index to be minimized 


K » Arbitrary constant associated with the 
penalty function technique 


= K for P 1 and P 2 , respectively 
L = Integrand of performance index 
m = Mass of movable object 
M = Mass of main body 

p = n column vector of influence functions 


P = Auxiliary state variable due to penalty 
function 



r = 


P for the two inequality constraints on 

x and x . , respectively 
max min r 

Position vector from center of mass of 
main body to the point mass 


r “■ Position vector from center of mass of 
main body to the center of mass of the 
system 


R q a Position vector from inertial origin to the 
u center of mass of the main body 

§ = Position vector from inertial origin to the 

center of mass of the system 


R = (n x q) matrix of influence functions 
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NOMENCLATURE (Continued) 


§ = First moment of mass of the system 

t = Time 

tg = Initial time 
tjr - Final time 

T = Rotational kinetic energy of the system with 
respect to its own center of mass 

u = m column vector of control variables 
u^ = Control variable equal to x 

W = (m x m) position definite matrix required in 
the first-order gradient method 

x 9 y,z = Coordinates corresponding to i 9 j 9 and k, 
respectively 

x s x « ” Maximum and minimum permitted mass positions 

I d ' x liuL " J on the x a :is 9 respectively 

x x - Actual values used for Xjnax and x m j_ n , 
n respectively j during computation 

x g = n column vector of state variables 

x s • = State variable 

X 7 Y,Z = Coordinate axes fixed at the center of mass 
of the main body 

3 “ State variable equal to x 

? = External moment 

e = Constant required in the first-order gradient 
method 

p = Equivalent mass equal to mM/ (M + m) 

v = q column vector required in the first-order 
gradient method 
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NOMENCLATURE (Continued) 


= Terminal constraint 


^1^2 = ^ for P 1 and P 2 5 respectively 

m = Angular velocity of spacecraft 

dj ,ta ,01 = Components of m along i, 5 > and ^5 

x y z respectively 

# ( 

til 1 ^ to , o)q = Components of u> along the maximum, inter- 
z mediate, and minimum inertia axes of the 

main body, respectively 

ui 9 , 01 « = Maximum values of oi 9 and to„ desired 

max d max 
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CHAPTER I 
INTRODUCTION 


General Introduction 

Future manned spacecraft might be subjected to 
tumbling. Such a situation may result from collision 
with another vehicle, thruster malfunction, propellant 
tank rupture, or escaping atmosphere. For example, the 
Salyut 2 space station orbited by the Soviet Union on 
April 3, 1973 went into a tumbling mode which is believed 
to be a result of an explosion or a wildly firing 
thruster. A study by Kaplan on tumbling causes showed 
that a future manned space station configuration may be 
subjected to a tumbling state with angular velocities 
up to two RPM if collision occurs with a space shuttle. 
Escaping atmosphere will yield about the same state. 

Tank rupture may result in angular velocities of approxi- 
mately one-half RPM. This tumbling should be immediately 
alleviated for crew safety and minimization of damage to 
the vehicle. Specifically, the crew would be subjected 
to oscillating accelerations . Hard docking by a rescue 
vehicle would be very difficult and would require a 
large fuel expenditure. In addition, it could not be 
implemented immediately. This time constraint also holds 
for other external detumbling methods such as fluid im- 
pingement. Hence, an internal control system is desirable. 
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Tumbling stabilization of space vehicles may be 
achieved by passive or active control devices. Passive 
systems are used as wobble dampers . They eliminate the 
wobbling motion of spacecraft by using mechanical or fluid 
dampers to dissipate energy until the minimum energy state 
is reached 5 their energy dissipation rates are low. This 
ultimately results in simple spin about the major principal 
axis. However, passive systems are designed for vehicles 
with high initial spin rates about the maximum inertia 
axis. Hence, they would not be appropriate for the future 
manned station mentioned above since it is normally in a 
non-spinning mode and a collision will result in the three 
angular velocity components of the vehicle being of the 
same order of magnitude. An active device such as a mass 
expulsion system may not be feasible since it requires 
long term, onboard storage of propellant. Momentum ex- 
change systems probably would saturate in large corrective 
maneuvers and may require continuous operation. An active 
control device that does not have the restrictions men- 
tioned above is an internal movable mass system. Movement 
of the mass will not affect the angular momentum vector, 
but it will affect the rotational kinetic energy. There- 
fore, internal mass motion can be used to decrease the 
rotational kinetic energy to a minimum which corresponds 
to the case of stable spin about the maximum inertia 
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axis. For the case of the previously mentioned future 
manned space station, there will still be need of de- 
spinning; but, docking by a rescue vehicle with subsequent 
despinning would then be a relatively simple operation. 

Also, escape by crewmen, if necessary, will be easier 
from a spinning rather than a tumbling vehicle. Movable 
mass control systems, however, have been investigated for 
vehicles which fall within the assumption of symmetry, or 
small transverse angular velocities relative to the spin 
velocity, or both. Stabilization of a vehicle like the 
manned space station cited above which is not symmetric 
and which may experience angular velocity components of 
the same order of magnitude requires further investiga- 

3 

tion. A recent study of this general problem by Edwards 
analyzed the rate of change of rotational kinetic energy 
in order to obtain stabilization; but, a long detumbling 
time and a large mass are required. Since it is important 
to have a mass as low as possible in space operations and, 
specifically for the distressed manned space station, to 
detumble as fast as possible, an optimal solution needs to 
be obtained. 

Statement of the Problem 

The objective of the research work presented here is 
to develop an optimal time and mass technique for obtain- 
ing the time history of internal control mass motion 
along a linear track in a tumbling space vehicle to achieve 
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simple spin. No assumptions of symmetry or small trans- 
verse angular velocities relative to the spin velocity 
will be made. Spec! fically , the manned space station 
mentioned in the previous section will be used as the 
test case. A first-order gradient optimization technique 
will be used to obtain motions of the internal mass that 
will result in stabilization about the maximum inertia 
axis. A penalty function method will be used to limit 
the extreme positions of the movable mass. Effects on 
the solution due to changes of various parameters will 
also be investigated. These parameters , all referring to 
the movable mass, are: mass, length of the linear track, 

positions along the two coordinates perpendicular to the 
linear motion, position along the axis of mass motion of 
the center of the track, and the direction of the track. 

A study of these changes will yield guidelines for obtain- 
ing maximum effectiveness from a movable mass control 
system. 

Summary of Work 

The differential equations of motion for a spacecraft 
with one internal movable mass permitted to move along a 
linear track were written with respect to an arbitrary 
orthogonal coordinate system fixed at the center of mass 
of the main body, which is the spacecraft without the 
control mass. A first-order gradient method, minimizing 

... 




the magnitudes of oscillations of angular velocity com- 
ponents along the intermediate and minimum inertia axes , 
was used to obtain control mass motions which yield 
simple spin about the maximum inertia axis. An IBM 370/165 
computer was used to obtain the results . This quantitative 
analysis, along with a qualitative examination of the 
differential equations of motion, permitted evaluation of 
various parameters in order to determine values which will 
result in minimum time and mass detumbling to simple spin. 
For minimum time, the mass, length of the linear track, 
and positions of the mass on the two coordinates perpen- 
dicular to the linear motion should have magnitudes as 
large as possible. Also, the mass should oscillate, 
about the zero point, on an axis parallel to the maximum 
inertia axis. A minimum mass solution is obtained by 
fixing the time at the largest feasible value. Compared 
to Edwards , this optimal technique permits detumbling in 
about one-fourth the time. Since stabilization may re- 
quire hours, this reduction in time is very significant. 

In regard to minimum mass, the optimization permits the 
use of a much smaller mass for detumbling in the same 
time. This mass reduction is quite substantial since 
very large masses are required. 




A, 


CHAPTER II 


PREVIOUS INVESTIGATIONS 

Differential equations governing the angular motion 

of a space vehicle with moving internal parts have been 

4 

obtained and discussed in several papers. Roberson de- 
rives these equations relative to the composite center of 
mass of the system. Since the reference point is the 
composite center of mass , the inertia dyadic is a function 
of time. The equations permit relative translational and 
rotational motion within the spacecraft frame. Effects 
on the vehicle due to known motions of the parts are then 
examined. Grubin 5 also obtains and discusses differen- 
tial equations governing the motion of a space vehicle 
with moving internal parts. However, his equations are. 
referenced to the center of mass of the vehicle without 
moving parts. Therefore, the inertia dyadic of the main 
body is constant. He examines several two dimensional 
cases of mass translation in a vehicle; but his examples, 
like those of Roberson, deal with effects on the vehicle 
due to known motions of internal parts . 

A number of papers in the literature deal with active 
control over the motion of internal parts in order to 
control the angular motion of a space vehicle. Since the 
differential equations are highly nonlinear, simplifying 
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assumptions are made in order to determine solutions . The 
vehicle is assumed to have small angular velocity com- 
ponents about two of its axes relative to the angular 
velocity along the axis about which final steady spin is 

desired, or is assumed to be symmetric, or both. 

n 

Kane and Scher analyzed the problem of active 
attitude control of a space vehicle with internal movable 
parts by considering its rotational kinetic energy. They 
noted that internal mass motion will not change the 
angular momentum vector since the net moment about the 
center of mass of the system is zero. However, internal 
mass motion will change the rotational kinetic energy of 
the system. Since the angular momentum vector is constant, 
rotational kinetic energy will be a maximum or a minimum if 
the rotation is about the minimum or maximum inertia axis 
of the vehicle, respectively. Thus, a tumbling space 
vehicle may be stabilized about the maximum or minimum 
inertia axis by internally moving a mass to either dissi- 
pate or add kinetic energy. If the space vehicle is sym- 
metric, they further noted that the kinetic energy of 
the vehicle may be used as a guide to determine the motion 
of an internal mass in order to have simple spinning motion 
or a combination of precession and spin in which the angle 
between the inertially fixed angular momentum vector and 
the axis of symmetry takes on any preassigned value. 
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Specifically, the kinetic energy is written in terms of 
this angle and the maximum and minimum kinetic energies . 

For a specified angular momentum vector, the initial and 
final kinetic energies may now be written by specifying 
the initial and desired final angle. A trial and error 
procedure is then used in order to find the motion of an 
internal mass that will result in a change of the initial 
kinetic energy toward the desired final value. This pro- 
cedure was applied to a solid uniform right-circular 

cylinder with a movable mass attached by a light rod. 
o 

Hopper investigated the use of internal mass motion 
in order to decrease the precession angle of spacecraft 
spin stabilized about their minimum inertia axis . Energy 
dissipating mechanisms are excited during precessional 
motion of this type of vehicle. This causes an additional 
increase in the precession angle since spin about the 
minimum inertia axis is one of maximum kinetic energy. In 
order to overcome this effect and have spin about the 
minimum inertia axis , energy must be supplied to the sys- 
tem. Two active devices are presented: one is a rotary 

device and the other is a linear oscillator. Both are 
examined for use in an axisymmetric spacecraft. The rotary 
device consists of a mass attached to and able to rotate 
about the minimum inertia axis . By keeping the rotor at 
some constant offset angle relative to its position due to 


... 
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centrifugal force resulting from precession of the 
vehicle, positive work can be done and thereby increase 
the kinetic energy of the system. The linear device 
consists of a small mass undergoing forced oscillations 
along a linear path fixed within the vehicle and per- 
pendicular to the spin axis. By oscillating the mass at 
the proper frequency and by proper control of phase, 
the driving motor will cause the kinetic energy of the 
spacecraft to increase . 

Childs investigated the problem of altitude stabili- 
zation of artificial-g space stations by a movable mass 
control system. The space station has the physical appear- 
ance of two rigid bodies connected by a long, slender tube 
and is spinning about its major principal axis. Movement 
by the crew may cause wobbling of the space station; that 
is, the angular momentum does not coincide with the maximum 
inertia axis. In order to damp this wobbling motion, the 
author placed a movable mass to the side of one of the two 
end pods. The mass was permitted to move inside of a tube 
which was parallel to the maximum inertia axis. By assum- 
ing small transverse angular velocities relative to the 
spin velocity, Childs was able to linearize the differen- 
tial equations of motion. This permitted the formulation 
of a control law to govern the motion of the movable mass 
in order to damp the transverse angular rates. However, 
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there was still some angular motion about one of the trans- 
verse axes after the application of the control law; the 
control system does not damp both transverse angular 
rates to zero. 

Lorell and Lange^ analyzed the use of internal 
moving masses to control the spin axis of a spinning 
satellite. Many satellites require control over the spin 
axis to an accuracy on the order of seconds of arc; two 
examples are communications satellites aiming high gain, 
narrow beam-width antennas and weather satellites scanning 
the surface of the earth for pictures. The control system 
has to take care of sensor-vehicle misalignments , motion 
of the principal axes of inertia, and body fixed disturb- 
ing torques . By assuming that the satellite is spinning 
about its axis of inertial symmetry, that the transverse 
angular rates are small relative to the spin velocity, 
and a specific geometry for four movable masses, the 
authors were able to simplify the differential equations 
of motion and use a linear control law. 

As stated in the introduction, Edwards performed an 
independent investigation of the general problem of a 
vehicle with arbitrary angular velocity components and 
arbitrary principal moments of inertia concurrently with 
this thesis. He formulated a control law for the force 
on a movable mass, in the direction of motion of the 
mass, that will reduce the rotational kinetic energy of a 

* ' . .(W _ [ x ^ ^ . . 
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tumbling vehicle and result in simple spin about the 
major principal axis. A formulation of a control law 
was made feasible by restricting the mass motion to lie 
along an axis parallel to the axis of maximum inertia 
and to pass through the zero point of that axis. The 
extreme mass positions permitted, on both sides of the 
zero position, are set in the control law; but, the 
formulated control law does not permit the mass to fully 
utilize the track available. Initially, the mass does 
go to an extreme position, but subsequent position peaks 
fall short of the extremes. Whether the extreme position 
will occur on the positive or negative side of the axis 
of mass motion is dependent on the initial conditions. 

Since these initial values are not known prior to instal- 
lation of a control system on board a space vehicle, the 
tube for mass motion must extend on both sides of the 
zero position a distance equivalent to the extreme posi- 
tion permitted. The detumbling times and masses obtained 
by Edwards are very large. For the manned space station 
mentioned previously, over three hours are needed to obtain 
stabilization using a movable object whose mass is 0.5% of 
the space station mass . An optimal time and mass analysis 
of the general problem is needed in order to reduce the 
mass and time required. 
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CHAPTER III 
ANALYTICAL STUDY 


Purpose and Procedure 

The purpose of this part of the investigation is to 
obtain the differential equations of angular motion for a 
space vehicle with a small movable internal mass and 
present an approach to obtaining mass motions that result 
in simple spin about the maximum inertia axis. Specific- 
ally, the rates of change of three orthogonal components 
of the angular velocity vector of the spacecraft will be 
expressed in terms of these angular velocities and the 
motion of the movable mass. For these equations, the 
center of mass of the main body will be used as the 
reference point for an arbitrary body fixed coordinate 
frame since this will allow the moments and products of 
inertia of the main body to be constant. Also, expressions 
will be derived for angular momentum, rotational kinetic 
energy, and rate of change of this energy. There will be 
no assumptions made of symmetry or of small angular 
velocities about two of the axes relative to a third; 
the space vehicle will be assumed to have three separate 
principal moments of inertia and an arbitrary angular 
velocity vector such that its three principal axis 
components may be of the same order of magnitude. 
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These equations are highly nonlinear and it is not 
obvious what the mass motion should be in order to detumble 
a space vehicle to simple spin about its maximum inertia 
axis. Whatever the motion should be, it must have reason- 
able values for mass displacement, velocity, and accelera- 
tion relative to the vehicle; that is, the mass should 
stay within the maximum dimensions of the space vehicle 
and should not have velocities and accelerations larger 
than can be supplied by a driving motor. Considering the 
above and further noting that a tumbling manned vehicle 
should be detumbled as quickly as possible, and possibly 
an unmanned vehicle from the standpoint of preserving 
structural integrity, the use of an optimization technique 
seems to be a feasible approach. Specifically, a first- 
order gradient optimization technique will b.i used since 
there are no previous solutions on which to base an 
initial guess of the control variable; this optimization 
technique does not require the initial guess to be close 
to the optimal values . Neighboring extremal and quasi- 
linearisation methods require good initial estimates of 
various parameters. Also, a penalty function method 
will be used to limit the extreme positions of mass 
motion. 
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Equations of Motion 

The angular momentum equation for a rigid body with 
an arbitrary origin of coordinates is 5 ’ 6 

? = ft + t x a Cl) 

where ? is the external moment , § is the angular momentum 
of the system, § is the first moment of md'5 of the 
system, and a is the inertial acceleration of the origin 
of coordinates . The desired equations of motion for a 
spacecraft with one small movable mass may be obtained from 
Equation Cl) by fixing the origin of coordinates at the 
center of mass of the main body, which is the spacecraft 
without the movable mass. The geometry of the system is 
shown in Figure 1 where x, y, z is a coordinate system, 
with 1 , ], k unit vectors, fixed in the main body, whose 
origin is at the center of mass of the main body. The 
angular momentum, 8, may be separated into two parts: 
the angular momentum of the main body relative to its own 
center of mass, fig, and the angular momentum of the mov- 
able control mass with respect to the center of mass of 
the main body, 8 . The angular momentum 8g may be 
expressed as 

g B = H Bx 1 + H By ! + H Bz 5< (2) 
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where 


Hr = 

D x 

I (0 - 

I <D - 

X „io„ , 

(3) 

X X 

xy y 

xz z 3 


% = 

- X to„ 

+ I m , - 


(4) 

xy x 

y y 

yz z 5 


H B Z = 

-I to 

-I to 

t I to . 

(5) 

xz x 

yz y 

z z 



The inertial time derivative of is 

+ “ x 

where [tj is the time derivative with respect to the 

£J 

body fixed x, y 5 z coordinate system of which is given 
by Equations (2) through (5) and to is the angular velocity 
of the spacecraft which can be expressed as 

to = to x i + oj y 3 + to^k. (7) 

The angular momentum consists of the angular momentum 
of the mass about its own center of mass and the angular 
momentum of the mass moving about the center of mass of 
the main body: 

§=$ +mrxr. (8) 

m m about its own 
center of mass 
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Considering the movable mass, m, as a point mass, the 
first term on the right hand side of Equation (8) can be 
equated to zero. The inertial time derivative of is 
then 

4 - . 

L = mr x r. (9) 

m 

The acceleration of the mass with respect to the center of 

11 

mass of the main body, r, may be expressed as 


r 


[r] tmxwxr + mxr , + 


2m x [r] 


where 


-*• 

r = 


** . “S’ 

xi + yj 


z£. 


Cr3 = xi + yj + z£, 

Cr] = xi + + zk. 


and, with a) given by Equation (7), 
m = (D x i + to y 3 + m z k. 

The first moment of mass of the system, is 


3 = 


mr . 


(10) 


( 11 ) 

( 12 ) 

(13) 


(14) 


( 15 ) 


The main body does not contribute since the reference 
origin is its own center of mass. Using Figure 1, the 
acceleration of the origin, a, may be expressed as 

- l t 
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where 



t 

M + m 


and, from the definition of center of mass. 


M 



mr 

M + m‘ 


(16) 


(17) 


(18) 


? is the resultant of external forces and M is the mass 
of the main body. Considering the spacecraft to be in a 
circular orbit with zero jet thrusting, F may be set 
equal to zero. Using Equations (16) through (18), a can 
now be written as 


r. S = - 


m 

M + m 


Noting that 



(19) 


(20) 


and placing Equations (9), (15), and (19) into Equation 
(1) gives 


? = a B 


+ 


mM 

M + m 


r x r. 


( 21 ) 




V 


u it. . 


& 




19 


For- zero external moment and noting that the force on the 
control mass is 


f = m(t + r) (22) 

m o 

which, upon using Equation (19) , becomes 

? = yr (23) 

m 


where 


U 


mM 

M + m 9 


(24) 


Equation (21) may be written as 

H n • = -r x ? . (25) 

B m 

Thus, as Equation (25) shows, the force on the control 
mass may be considered as causing a moment to act on the 
main body. The mass will be permitted to move along a 
linear track parallel to the x axis. Placing Equations 
(2) through (7), (10) through (14), and (23) into Equa- 
tion (25) yields three scalar equations which, when 
solved simultaneously, give the differential equations 
of motion in the desired form: 
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01 = d) (ul , W , Ui , X, X, x) , 

x x x 5 y 5 z 5 s 5 5 


( 26 ) 


<i y = d) (u> x9 w y , o) z , x, x, x), 


( 27 ) 


d) = (I) (d) , 0) , 0) , X, X, x). 

z z x’ y 5 z ’ ’ 5 


( 28 ) 


The full equations for ui^, w , and are given in 
Appendix A. Since the motion is along an axis parallel 
to the x axis , 


♦ ** • *• __ ft 

y“y-z = z — 0 


( 29 ) 


and the y, z positions of the mass may be arbitrarily 
fixed. Other constants that have to be specified are 
moments and products of inertia, and masses of the main 
body and the movable object. The control mass was per- 
mitted to have motion parallel to just one axis since then 
only one control variable will need to be specified; this 
fact will become important in the following section which 
deals with the optimization technique. However, the direc- 
tion of the x axis, and, therefore, the direction of mass 
motion, relative to the spacecraft may be changed arbi- 
trarily by appropriately changing the moments and products 
of inertia in the equations of motion given in Appendix A. 

The total angular momentum vector with respect to the 

center of mass of the system, 3 . must remain constant 

cm * 
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during control mass motion. 


that 


11 


This can be seen by noting 


r 



(30) 


and, since there is no external moment during mass motion, 

H is constant. An expression for this total angular 

momentum vector can be obtained by dividing it into two 

parts : that due to the rotation of the main body and that 

due to the motions of the centers of mass of the main body 

and the movable object about the center of mass of the 

system. The former is simply I * m where I is the inertia 

-y 4 - 

dyadic. The latter can be expressed as pr x r by con- 
sidering the two-body problem composed of the main body 
and the movable mass as an equivalent one-body problem 

which is the equivalent mass p moving at a distance r 

11 

from the center of mass of the system. Thus, we have 

S=I*S+prxr (31) 

cm 


where 


4 - 

r 



+ w 


X r. 


(32) 


With this total angular momentum vector constant , the 
rotational kinetic energy, T, can assume the following 
values : 
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(33) 


where I and I . are the spacecraft ' s maximum and min- 
max mm r 

imum moments of inertia. Specifically , using the same 

analysis as was used to obtain H , we can write the 

cm 

following expression for kinetic energy relative to the 
center of mass of the system during mass motion: 

T = • I • w + ^yr * r. (34) 

To have simple spin about the maximum inertia axis, the 
tumbling vehicle’s rotational kinetic energy, which ini- 
tially is some constant value, must be decreased to the 

2 

value associated with this simple spin, H /2I Q . It 

r r 7 cm max 

should be noted that simple spin about the major principal 
axis of the spacecraft can essentially be considered, 
in this investigation which uses a small mass, as simple 
spin about the major principal axis of the main body. 

If the control mass is large and far from the center 
of mass of the main body, the orientation relative to 
the spacecraft of the maximum inertia axis of the main 
body may be quite different from that of the spacecraft. 

In that case, the maximum inertia axis of the spacecraft 
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should be used when referring to simple spin about the 
maximum inertia axis . Since the change in the kinetic 
energy of the system composed of a rigid main body and a 
movable mass requires work to be done and the only source 
of work is the force ¥ acting on the control mass m 
which moves a distance dr, we have 

d(Work) = ¥ * dr. (35) 

m 

This equation can also be obtained by considering ¥ 
which is given by Equation (23) as acting on the equivalent 
mass m and causing a displacement dr. For mass motion 
along a linear track parallel to the x axis , the right 
hand side of Equation (35) may be written as f mx dx. 
Therefore, we can write the following expression for 
the rate of change of kinetic energy which is equal to 
the rate at which work is being done : 

T = f mx *. (36) 

It should be noted that T given by Equation (36) determines 
the power that will be required. 

Optimal Control 

A first-order gradient algorithm for the following 

12 13 

problem is available in the literature: s with u and 

x defined as the column vectors of the control and state 
s 

variables, respectively, u^(t), ... u m (t) must be found 
in order to minimize 


24 



J = / L[Xg^(t), . .. Xg n <t) S 

t o 


U l (t) ’ ••• U IE (t) 


where 


t ]dt 
(37) 


“ F j/ X si j * * * x s n > 5 * • • U jn 5 i-1, * * • n 

(38) 

with x g (t Q ) , t QJ t f specified, and terminal equality con- 
straints on q of the x s . variables, each represented by 

JL 

l|j[x s . (t f )] = Xg. <t_) - Xg. = 0 (39) 

1 r 1 r ^final 

with the q desired terminal values, x s . , specified. 

1 f inal 

Inequality constraints of the form 

W u, t) < Fp max (40) 

may be handled by using a penalty function technique 
which converts inequality constraints to terminal con- 
straints . Specifically, an auxiliary state variable P 
is defined as 


p = { KCF P Cx s ’ u > “ ^fixed 3 V F P - Fp fixed (41) 
° ! F P < Fp fixed 




U _ 


and then the state variable P is forced to approach a 
value of zero by specifying the following form of Equa- 
tion (39): 

ip = P(t f ) = 0. (42) 

In Equation (41), K is an arbitrary constant, and Fp^ ^ 
is chosen to be smaller in magnitude than Fp max since, 
in the first-order gradient method, P needs to exist in 
order to be controlled. If ^ is set equal to 

Fp_ 3 then some violation of the inequality constraint 
will have to be accepted. 

The problem specified in the previous section will 
now be put into the form required for the application of 
the first-order gradient method. The dif ferential equa- 
tions of motion given by Equations (26) tl rough (28) can 
be put into the form of Equation (38) by the following 
substitutions : 

& = g (43) 

♦ 

and, since x is equal to g, 

$ = u 1 . (44) 

The state variables, x s . with i = 1, ... 5, are m , u . 

x x y 

m , x, and g, respectively. The one control variable is 
z 

u^j more control variables would be needed if the mass 


was not restricted to move parallel to one axis . The 
equations of motion can now be written as follows : 


^x = 

“x <l V 

V 

to , 
z 5 

X, 

e, 

s 

(45) 

*y = 

V“x’ 

“y> 

z 5 

X, 

3, 

u l ) J 

(46) 

0) = 
z 

“z (( V 

V 

0) , 

X, 

3, 

Ui), 

(47) 

X = 

e, 






(48) 

• 

B = 

u r 






(49) 


In order to have simple spin about the maximum inertia 
axis, co 2 and which are angular velocity components 
along the intermediate and minimum inertia axes need to be 
minimized* Thus, modeling after regulator problems, the 
performance index J given by Equation (37) will be ex- 
pressed as 



w, 


to, 


•)dt 


to. 


w. 


‘max 


max 


(50) 


where and m 3 are the maximum magnitudes that are 

desired. Ideally, these values should be zero to have 
pure, simple spin about the major principal axis. However 
in practice, these maximum values will be set at some very 
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small magnitudes. The variables and , of course, 
need to be expressed in terms of the state variables. 

If, for example the x, y, and z axes are aligned with the 
maximum, intermediate, and minimum inertia axes, respec- 
tively, then w 3 are equal to and a> z . 

The x position of the control mass is limited as follows: 


x . < X < X 

mm - max 


C 51 ) 


with the values x . and x arbitrarily set. In order 

mm max 

to apply the penalty function technique. Equation (51) 
will be expressed as 


x < x 


max 


and 


-x < 


x . . 

mm 


(52) 


(S3) 


Thus, there will be two auxiliary state variables, 
and P ^ , associated with Equations (52) and (53), respec- 
tively; state variables x S g and x s ^ will refer to and 
P^. Specifically, we have the following additional 
equations of motion: 


p x -- { 


V x - V * x - *h 


(54) 


0 , x < x T 
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and 


P 


2 


K (~x + x 1 ) 2 , 
{ 2 1 
0 5 -x < -X-j 


-X > -X x 


(55) 


where x, is less than x and x n is greater than x . . 
h max ± mm 

Using Equation (42), we can now place terminal constraints 
on the two auxiliary state variables as follows : 

i|>l = Pj^Ctf) = 0 (56) 

and 

i|» 2 = P 2 <t f ) = °* (57) 

The steps of the first-order gradient optimization 
technique for the problem just specified can be written 
as: 12 ’ 13 

1. Estimate u-^Ct) and integrate Equations (45) 
through (49), (54), and (55) forward with the known 
initial conditions x g (tg). Store x s (t), u^(t), tJv 1 , 
and if) . 

2. By backward integration of the following influence 
function equations determine and store the n column 
vector p(t) and the (n x q) matrix R(t) where n equals 
seven and q equals two: 
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P = 


-< 


I^-) T P - & ) T 


(58) 


and 


R 


,3f > T r 

-‘aST 5 R ’ 


(59) 


with 

pCt f ) = 0 


(60) 


and 


9 


1, 


R i3 C V = 


0 

1 


s 
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i = 
i = 

i - 

i = 
i = 


1, . . . 5 and j 
6 and j = 1 

6 and j = 2 

7 and j = 1 

7 and j = 2 


1 , 2 


(61) 


where 1=6 and 7 correspond to the auxiliary state vari- 
ables and respectively. 

3. Evaluate the following integrals: 


R dt 


(q x q) matrix, 


(62) 


& = i‘ CpT + W " 1( % )T R dt 


q row vector. 




A, 


(63) 
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/ c p T * |i_) w _1 [(|£-) T p + 


3u. 


3u. 


3u. 


3u- 


scalar , 


(64) 


where W is just a (lxl) positive-definite matrix since 
there is only one control variable. 

4. Select a 6 i|j, given by 


StJ> 

fill; = l 1 1 9 (65) 

64, 2 


which will bring i|>^ and ip 2 , given by Equations (56) and 
(57), closer to zero on the next iteration. Specifically, 
choose 

fii|f^ =■ eP^(t^) (66) 

and 

2 = eP^t^) (67) 

where 

0 < e < 1. (68) 


Then, determine the q column vector v from 


v 


< 6 * + v 


(69) 
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5 . Repeat steps 1 through 4 with the following 
improved u^Ct): 


u l (t) new 


= u l (t) cld + 6u l Ct) 


( 70 ) 


where 


<$Ui (t ) 




3u. 


+ CpCt) + R(t)v 


5 f 
9u. 


-]} 


(71) 


T -1 

Terminate when and Ijj - equal zero 

to the desired degree of accuracy. These steps have been 
speciliazed to the problem consisting of seven state 
variables , with terminal constraints on the sixth and 
seventh state variables, and one control variable. 

Numerical Solution 

The five steps of the first-order gradient optimiza- 
tion technique were programmed on an IBM 370 /IB 5 computer. 
The computer program, given in Appendix B, consists of a 
main program with three subroutines* In the main program, 
the necessary variables are specified. The variables 
needed for the specification of a case which is to be 
studied are listed in the beginning of Appendix B, their 
computer language names are also specified. The sub- 
routines carry out the integrations -and changes in the 
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control variable as prescribed in the five steps . The 
differential equations are solved by using a fourth-order 
Runge-Kutta algorithm. Integrals are handled by using 
a library program consisting of an extended five-point 
Newton-Cotes quadrature formula. 

Implementation and Nature of the Optimal Control 

The optimal solution obtained is really a local 
optimal. Standard numerical optimization techniques like 
the first-order gradient do not necessarily yield the 
absolute minimum. Initial guess of the solution will 
determine which local optimal is obtained ; this optimal 
may be the absolute optimal. Of course, for a specific 
problem there may be only one minimum and, therefore, the 
solution Is the absolute minimum. An examination of 
numerical results and a comparison to non-optimal solu- 
tions will give an indication of the nature of the solu- 
tion obtained. Results of this investigation are pre- 
sented and compared to those of Edwards in the next 
chapter . 

The optimal control may be implemented quite easily. 
The first-order gradient method yields a time history of 
the mass motion. Thus, the control system need just 
monitor and change the position of the mass. The position 
may be monitored by a simple mechanical device. It should 
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be noted that any means may be used to move the mass. In 
actual space operations 5 three orthogonal angular veloci- 
ties of a tumbling vehicle will be sensed by rate gyros 
and extrapolated to a time a few minutes in the future 
using Euler moment equations. These future angular 
velocities will then be used as the initial conditions 
for the optimization equations . The optimal control will 
be numerically obtained and initiated at the chosen 
future time. 
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CHAPTER IV 
RESULTS 

Several parameters need to be specified before 
computer simulations may be run for minimum time stabiliza- 
tion. As will become evident, selection of some of these 
parameters will depend upon the specific satellite to 
which the movable mass control system is applied; that is, 
the dimensions of the satellite and the amount of time 
that can be permitted for detumbling. The choices for 
minimum time detumbling of the remaining parameters will 
also become evident; but, these parameters will not be 
dependent on the type of satellite to be controlled. 
Specifically, referring to the differential equations 
of motion in Appendix A and noting that the mass moves 
along an arbitrary x direction, these parameters are as 
follows: mass of the movable object, length of the 

linear track, y and z positions of the mass, point about 
which the mass oscillates , and direction of the x axis 
relative to the spacecraft. By examining Equation (25) 
and thinking in terms of moments applied about each axis , 
a qualitative preliminary analysis may be made as to the 
effect of these parameters on the time needed to detumble. 
Increasing the mass of the object will increase the force; 
thereby increasing the moment and permitting a decrease in 
the detumbling time. Increasing the length of the linear 
track or the y and z magnitudes will increase the moment 

► * . ... u A, _ . _ 
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arm, which should tend to increase the moment and result 
in a lower minimum time. It should be noted that changes 
in these parameters will also effect the force; but, the 
overall effect on the moment will probably be due to the 
change in the moment arm mentioned previously since the 
force consists of the relative difference of various 
terms. By increasing the x value of the point about 
which the mass oscillates- the moment arm is again in- 
creased. However, by permitting the mass to move further 
on one side of the zero x position than on the other, 
there may arise difficulties due to a larger moment in 
one direction than in the opposite. Changing the direction 
of the x axis relative to the spacecraft will affect the 
moment arms of the force components producing moments 
about the intermediate and minimum inertia axes. If the 
x axis is parallel with the maximum inertia axis , maximum 
control over the moment arms of the moments about the 
intermediate and minimum inertia axes will be available; 
this can be seen by noting that now the x axis is per- 
pendicular to the intermediate and minimum inertia axes . 
This would seem to indicate that the orientation of the 
linear track should be parallel to the final spin axis , 
which is the major principal axis. However, in this case 
as in the other cases that involved changes in the moment 
arm, there are also changes in the force itself which are 
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difficult to specify. The above qualitative analysis 
is not sufficient in itself to arrive at minimum time 
values for all the parameters being discussed. A 
quantitative analysis must be made. 

The movable mass control system will be applied to 
a manned space station configuration which NASA is con- 
sidering for the 1980’s. This configuration is shown in 
Figure 2 with pertinent data being given in Table 1. 

The optimization procedure will permit the fastest de- 
tumbling possible with the movable mass. Based on the 
previous qualitative discussion of various parameters in 
the differential equations of motion , we will initially 
fix the parameters to yield the best detumbling sequence; 
further cases will vary these parameters in order to 
quantitatively determine the minimum time solution. The 
mass of the movable object will be set at 499 kg, which 
is 0.5% the mass of the manned space station. It will be 
permitted to travel approximately ±3.7 m about the zero 
position on the x axis. This axis of motion will be 
parallel to and have the same sense as the maximum 
inertia axis. For conveniences the y and z axes will be 
chosen to be parallel to and have the same sense as the 
intermediate and minimum inertia axes, respectively. 
Choosing the y and z positions are large as possible 
within the limits of the space station, we have 5.55 m 
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Table 1. 


Manned Space Station Data 


2 9 15 


Mass 

Maximum Moment of Inertia 
Intermediate Moment of Inertia 
Minimum Moment of Inertia 


99792 kg 
; .74 x 10 6 kg-m 2 
6.28 x 10 6 kg-m 2 
5,15 x 10^ kg-m 2 


Transformation matrix, body fixed X, Y a 2 to principal 1, 
2, 3, where 1, 2, 3 are the maximum, intermediate, and 
minimum moments of inertia axes respectively 


0.458 

-0.889 

0.00676 

0.889 

0.458 

-0.00204 

-0.00128 

0.00695 

1, 0 









& 


39 


and -13.7 m. The initial angular velocity components 
along the 1-axis of maximum inertia, the 2-axis of inter- 
mediate inertia, and the 3 -axis of minimum inertia will 
be chosen as 0.103 rad/sec, -0.199 rad/sec, and 0.000286 
rad/sec. These values are based on a worst case tumbling 
analysis. They represent the highest tumbling mode of 
the manned space station and are due to a collision 
between it and a space shuttle vehicle. If uncontrolled, 
the manned space station with a fixed 499 kg internal 
control mass would continue to tumble with oscillating 
between 0.103 rad/sec and 0.192 rad/sec, tog oscillating 
between -0.199 rad/sec and 0.199 rad/sec, and u>g oscillat- 
ing between -0.118 rad/sec and 0.118 rad/sec. Flexibility 
effects, of course, will tend to decrease the rotational 
kinetic energy of the vehicle and alter the envelopes of 
oscillation. However, in the time periods involved in 
control, these limits of oscillation can be used as a 
reference for zero control mass motion. The effects of 
control by an internal movable mass system on the oscil- 
lations of to^, togs and are shown in Figure 3. The 
curves in this figure are the envelopes of oscillations 
of the principal axes angular velocity components. At 

2,845 sec, the limits, for the penalty function, of 

-9 

mass motion along the x axis were set at ±10 m m 
order to zero out the mass position, velocity, and 
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acceleration. By 2,893 sec these three variables were 
essentially zeroed out, having values of 0.0466 m, 

-0.00322 m/sec, and 0.00141 n/sec . After this time, the 
mass was kept fixed at the zero x position*, remained 
at 0.212 rad/sec, oscillated between -0.00152 rad/sec 
and 0.00142 rad/sec, and oscillated between -0.000837 
rad/sec and 0.000909 rad/sec* It should be noted that the 
percentage change from peak to peak of and is greater 
toward the end of the control time since the absolute de- 
crease remains fairly constant. After about 3,000 sec, 
then, the peak values of tUg and were reduced by more 
than 99% of their initial value. Of course, the mass con- 
trol system could have been left on to reduce the and 
oscillations even further. For the manned space station 
discussed in this investigation, these values are suffi- 
cient since the effect felt by the crewmen is essentially 
that of simple spin and docking by a rescue and despinning 
vehicle can be made as if with a simple spinning body. 
Figure 4 shows the motion of the intei'nai movable mass 
which results in detumbling of the manned space station. 

The mass position, x, at no time exceeds 3.7 m; after 
2,893 sec, it is essentially equal to zero. The mass 
oscillates between its maximum permitted limits , utilizing 
the full linear track available to it. The velocity of 
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the mass, x, oscillates between -0.647 m/sec and 0.654 
m/sec. The greatest mass acceleration, x, occurs during 
the zeroing out of the mass position, velocity, and ac- 
celeration: in this period it reaches values of -0.579 
m/sec and 0,465 m/sec . Also during the zeroing out, 
the force in the x direction, f m , acting on the mass 
reaches its largest magnitude, 283 N. These values for 
mass velocity and acceleration, and for force, are reason- 
able. It should further be noted that the mass velocity 
and force maximum magnitudes occur during energy dissipa- 
tion, -T , During energy dissipation, the force in the x 
direction and the mass velocity are opposite in direction; 
here the control system is actually restraining the mass. 

7 

As Kane and Scher have pointed out, this energy dissipa- 
tion could be used co provide useful power for the vehicle’s 

systems. Energy hss to be provided by the mass only when 

• • 

T is positive. Figure 5 shows the variation of T, power 

with time. The positive power is much less than the nega- 
tive. Much more energy is dissipated than added. Also, 
as this figure shews, the total energy that is to be 
supplied is quite reasonable. The maximum positive power 
is 48.4 watt, also a reasonable value; this value could 
have been reduced to a value similar to the other peaks 
by decreasing the simulation time increment. Time incre- 
ments will be discussed later. The -30G watt value at 
2,848 sec, during zeroing out, corresponds to energy 








dissipation; the mass is being restrained and the control 
system does not have to supply that energy. Figure 6 
shows the decrease of rotational kinetic energy from its 
initial , before control mass motion, value of 1.62 x 10 5 
joule to the final value for stable simple spin of 

c 

1.5 x 10 joule. At various points in this figure, 
the kinetic energy Increases slightly and then resumes 
its downward curve. These increases, of course, cor- 
respond to energy addition to the movable mass control 
system, the positive power points on Figure 5. Further- 
more, by superimposing Figure 4 for mass x position with 
Figure 5 for power during mass motion, it is seen that 
these energy addition points correspond to the points at 
which the mass direction of motion needs to be reversed 
in order not to exceed the extreme limits of motion that 
were previously set. Throughout the period of control 
mass motion, the angular momentum of the system relative 
to Its center of mass remains fixed at the value it had 

g 

before control mass motion was initiated, 1.45 x 10 
kg-m /sec. The angular momentum vector remains constant 
since there are no external moments on the space vehicle. 
Therefore, as is evident, an optimal movable internal 
mass control system can be used to reduce the arbitrary 
tumbling of a general space vehicle to simple stable spin 
about the maximum inertia axis . 


Kinetic Energy t^oule; 


1.62x10 


1. 59x10 


1.56x10 


1 . 53x10 



Figure 6. Decrease of Rotational Kinetic Energy Using 0.5% Mass 
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One valid comparison between the optimal control 
mass motion described in this investigation and the limited 
analysis of the force control law formulation would be 
to keep all parameters and initial conditions mentioned 
in the previous paragraph the same , including the per- 
mitted extreme positions on both the positive and negative 
sides of the mass motion axis; as stated in the previous 
investigation, the linear track required by the force 
control law method must extend, in the positive and nega- 
tive directions, a distance equivalent to the maximum 
movement permitted* Doing this, the force control law 
yielded a decrease of the envelope of oscillation to a 
magnitude of 0.00206 rad/sec at 11,050 sec, and the w 3 
envelope to a magnitude of 0.00124 at 11,005 sec. As was 
mentioned in the previous paragraph, the optimal analysis 
of this investigation decreased the magnitudes of the 
and 0)^ envelopes to 0.00152 rad/sec and 0.000909 rad/sec 
respectively by about 2,900 sec. Thus, the optimal 
analysis permits detumbling in approximately one-fourth of 
the time required by the force control law analysis. 

Since the one-fourth value means that the crewmen will be 
subjected to a tumbling state for less than an hour 
compared to over three hours, it is quite significant. 

The force control analysis required only about one watt 
of peak positive power. As was shown, the power for the 

L. ^ ... 
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optimal control is greater, but still within the limits 
of production in spacecraft. Another comparison between 
the optimal analysis and the formulated force control law 
approach can be made by restricting the mass motion 
extreme positions in the former method to the actual 
extreme positions of the latter, keeping all other ini- 
tial conditions and parameters the same 5 the values used 
in the case mentioned in the previous paragraph will 
again be used. The first 160 sec after commencing 
control mass motion were investigated. In this time 
interval, the force control law required mass position 
peaking of 0.631 m at 35 sec and -3.75 m at 155 sec; 

^2 and peak at 0.19 6 2 rad/sec and 0.1179 rad/sec. The 
next peaking of mass position occurs at 295 sec with a 
value of 2.23 m. The optimal mass control system was 
started at the 50 sec point of the force control sequence 
since, by that time, the mass positi.on had peaked at 
only 0,631 m. Using the values at the 50 sec point, the 
optimal method peaked the mass position first at the 
largest positive limit and then at the negative one. 

This motion resulted in an ( 1)2 peak of 0.1947 rad/sec 
and an a > 3 peak of 0.1171 rad/sec. Without any control 
mass motion, the vehicle would have experienced an ui 2 
peak of 0.1991 rad/sec and an to 3 peak of 0.1182 rad/sec. 
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Thus, the optimal technique yielded about a one and one- 
half times greater decrease in the value and slightly 
less than four times greater decrease in the value. 

The actual effect on the time to detumble will be greater 
since the angular velocities are lower for the start of 
their next cycles and their periods are decreasing more; 
this is in comparison to the values yielded by the force 
control law method. Therefore, the optimal control 
mass motion technique, in addition to not having the 
restriction on the direction of mass motion nor on the 
point about which the mass oscillates , yields simple spin 
in a considerably faster time than the force control 
law method. 

The effect of a change in the mass of the movable 
control object was investigated by using the case initially 
studied and only changing one parameter, the mass, from 
499 kg to 998 kg; this new value for the mass is one per- 
cent that of the manned space station. Specifically, the 
control mass will again be permitted to move 3.7 m about 
the aero position on the x axis which is parallel to and 
has the same sense as the maximum inertia axis, with the 
y and z positions fixed at 5.55 m and -13.7 m, respec- 
tively. Also, the y and z axes are parallel to and have 
the same sense as the intermediate and minimum inertia axes. 
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The initial values for w.. , and ui„ remain set at 0.103 

i. I o 

rad/sec, -0.199 rad/sec, and 0.000286 rad/sec. Figure 7 
shows the envelopes of oscillation of the principal axis 
angular velocity components. Comparing with Figure 3, 
it is seen that doubling the mass to one percent of the 
space station mass caused about a one-half decrease in 
the time to detumble to simple spin about the maximum 
inertia axis. After 1,485 sec, to 2 is oscillating 
between -0,00186 rad/sec and 0.00199 rad/sec, and 
between -0.00126 rad/sec and 0.000425 rad/sec. To reach 
similar oscillation ranges, the 0.5% mass required about 
3,000 sec. Comparing this one percent case to Edwards, 
it is again evident that only about one- fourth the time 
is required; his force control law reduced u 2 and u > 3 
to peak magnitudes of 0.00292 rad/sec and 0.0032 rad/sec 
at 5,680 sec and 5,450 sec, with mass set at one percent. 

The effect of a change in the length of the linear 
track is evident in every computer run. At first, the 
peaks of the m 2 and oscillations are lowered in magni- 
tude considerably by having the mass move far out on the 
x axis. Each successive iteration improves on the control 
variables in order to decrease the extreme positions of 
the mass to prescribed limits. As the extreme limits are 
decreased, the peaks of the ti > 2 anc ^ ^3 oscillations increase 
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in magnitude. The above examined case of a one percent 
control mass will be used to show this effect. The first 
100 sec of simulation are shown in Figure 8. After a 
few iterations , the mass position extremes are reduced 
to ±10 m. The resulting oi^ s Wg » and oscillations 
are plotted. After additional iterations, the mass posi- 
tion limits are -2.B m and 3.4m, and it is seen in this 
figure that the and oscillations have become 
larger in amplitude. Comparison to the no mass motion 
curves shows that the effect of extreme mass position 
change is quite substantial. It should be noted that, as 
the extreme limits of mass motion are decreased and cause 
an increase in the amplitudes of the and oscilla- 
tions, an increase in total time to detumble will occur. 

The changes in the oscillations of and due to 
a one-half decrease in the y and z positions of the control 
mass are examined by using the one percent mass case men- 
tioned above and appropriately changing y and z, A 998 
kg mass was permitted to move 3,7 m about the zero posi- 
tion on the x axis which is parallel to and has the same 
sense as the maximum inertia axis, with y and z fixed at 
5.55/2 m and -13.7/2 m, the y and z axes are parallel to 
and have the same sense as the intermediate and minimum 
inertia axes. The only differences between the case that 
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will be discussed now and the initially discussed 0.5% 
case are an increase in the mass of the movable object to 
one percent of the manned space station and a one-half 
decrease in the y and z positions which were 5.55 m and 
-13.7 m. Figures 9, 10, and 11 show the a) 2 , and io 3 
oscillations. For comparison, the o) 2 and oscillations 
from the previous one percent case, which had the same 
parameters except for full y and z values of 5.55 m and 
-13.7 m, are shown. The limits of the w 2 and oscilla- 
tions for no motion of a one percent mass are also noted. 

By comparing the various curves , it is seen that there 
is a definite increase in the magnitudes of the and to ^ 
oscillation peaks caused by lowering the y and 2 mass 
position magnitudes. The long term effect will be an in- 
crease in total detumble time due to the lower y and z 
magnitudes . 

The effect of a control mass oscillating about a 
non zero position is shown in Figure 12. The initial 
one percent mass case was again used as a basis and 
the only parameter change was letting the mass oscil- 
late about +10 m instead of the aero x position. Specific- 
ally, a 998 kg mass was permitted to move 3.7 m about the 
10 m position on the x axis which is parallel to and has 
the same sense as the maximum inertia axis, with y and 2 
fixed at 5.55 m and -13.7 m; the y and z axes are parallel 
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Figure 10. t02 Oscillations Using One Percent Mass with Various Values for 

y and z 







Figure 11. ( 1)3 Oscillations Using One Percent Mass with Various Values for 

y and z 
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to and have the same sense as the intermediate and minimum 
inertia axis. Initial values for m 2 , and Wg are 

still 0.103 rad/sec, -0.199 rad/sec, and 0.000286 rad/sec, 
respectively. For both the a) 2 and curves, it is evident 
that one side of the oscillation is decreased faster than 
the other. One side of the envelope crosses the time axis 
and, along with the other side, tends to a small, but 
finite, value. The co^ and Wg oscillations, then, will 
not be completely zeroed out*, but, they will be made quite 
small. Comparison to the case shown in Figure 7 which 
differs only in that the mass oscillation is about the 
zero x position instead of 10 m, shows that the higher 
x values initially permit a much faster decrease in the 
and u> 3 envelopes of oscillation. However, having the 
control mass move through the zero position on the x axis, 
results in the u> 2 and tOg oscillations tending to zero 
rather than a finite value. 

The movement of the control mass parallel to an axis 
other than that of maximum inertia was investigated by 
having the mass move parallel to the intermediate inertia 
axis . This choice of axis and that of the other parameters 
was made in order to let this case be as similar as possi- 
ble to the initial one percent mass case; this permits a 
more valid observation of the effect of the change in the 
direction of control mass motion relative to the main 




60 


vehicle. The problem arises since a change in the x 
direction of mass motion changes the y and z values of the 
mass and what they represent relative to the main vehicle. 
Therefore j for this case under consideration, a 998 kg 
mass was permitted to move 3.7 m about the zero x posi- 
tion. The x axis , along which the control mass moves, was 
placed parallel to and has the same sense as the inter- 
mediate inertia axis. The y and z axes were placed 
parallel to and have the same sense as the minimum and 
maximum inertia axes. The y and z positions of the mass 
were set at -13.7 m and 5.55 m. Initial values of w to^, 
and tog remain at 0.103 rad/sec, -0.199 rad/sec, and 
0.000286 rad/sec. Oscillation envelopes for to , w^, and 
are shown in Figure 13. The to^ envelope tends to zero. 
However, the to^ envelope is tending to a small, but finite, 
value. Comparing i,e case shown in Figure 7 it is 
seen that a cont. o' ma.st moving parallel to the maximum 
inertia axis permits , in .Addition to zero values for 
both the to^ and tog envelopes, faster detumbling. 

These results were obtained by running the comptuer 
program for simulation times of 100 sec. To bring the 
extreme mass positions to the permitted magnitudes , in 
this 100 sec simulation time period, required up to about 
35 iterations which used about 100 sec of IBM 370/165 
computer time. At the end of each 100 sec simulation 
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Figure 13. Envelopes of Angular Velocity Oscillations Using One Percent 
Mass Motion Parallel to the Intermediate Inertia ^xis 






x and x were used as 


run the end values for o)g , 

initial conditions for the next 100 sec simulation run. 
This time increment of 100 sec for simulation was chosen 
since larger time increments made it difficult to bring 
the extreme mass positions to within the specific limits. 

A few times even this time increment was too large ; that 
is, the mass would initially stay within the limits but 
would then move past, where about 35 iterations were the 
maximum permitted. Rather than extend the computer time 
and thereby increase the number of iterations , the accept- 
able initial part of the run was used. For zeroing out 
the mass position, velocity, and acceleration in the 0.5% 
case that was initially investigated, the simulation time 
increment was arbitrarily chosen at 50 sec. Also, the 

-9 

position limits, as stated previously, were set at ±10 m 
for the zeroing out simulation time. Normally, limits 
were set at ±2.5 m at the beginning of a case and then 
changed to about ±3.0 m or higher. As stated previously 
in the analytical study on the optimisation technique, 
the limits of mass position should be set lower than what 
is desired since the penalty function comes in when there 
is a violation of the set limits . The iteration that was 
chosen during each 100 sec simulation time run had the 
lowest peak magnitudes for the tOg and (0g oscillations for 
mass positions within the ±3.7 m prescribed extreme 
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limits. There was no need to place constraints in the 
optimization technique on mass velocity and acceleration 
since, as discussed in the initial 0.5% mass case, these 
variables did not reach excessive magnitudes. The values 
of other constants, associated with the optimization method, 
which were discussed in the analytical study are given in 
the main program of the computer program listed in 
Appendix B. The time steps in the integrations were set 
at 5.0 sec of simulation time. 
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CHAPTER V 

CONCLUSIONS AND RECOMMENDATIONS 

An optimal movable mass control system has been 
applied to a tumbling spacecraft in order to obtain simple 
spin about the major principal axis. The results indicate 
that the largest possible magnitudes for the internal mass, 
length of the linear tracks and positions of the mass on 
the y and z axis will yield the fastest detumbling times. 
The choice of these values depends upon size and mass of 
the spacecraft. Results also indicate that the mass 
should oscillate s about a zero point, on a line parallel 
to the maximum inertia axis. These results were based on 
worst case initial conditions for to., , and . 

Tumbling situations that might be encountered in actual 
space operations will usually be less severe and, there- 
fore, will probably require less time to reach simple 
spin. Also, these results were based on one vehicle, 
the modular space station. However, since this vehicle 
was asymmetric, the optimization technique will apply 
to any type of spacecraft. The results of various para- 
meter changes, furthermore, were based on motions of a 
one percent mass and the resultant effects on the peaks of 
and cd 3 i no further comments were made about T, 
f , x, and x. A one percent mass was used to show the 
effects of changes of various parameters since the large 

L . & ... , . 




65 


mass made the effects readily apparent and showed them in 
a faster time, compared to smaller masses which may be 
more feasible for the space station due to its considerable 
mass. The graphs of 5 , and were used since the 

objective is to reduce the peaks of o) 2 and io 3 , with 
tending to one value. Other variables were not discussed 
since their behavior and magnitudes were comparable to 
the 0.5% mass case which was examined in detail. This 
0.5% mass case showed that the velocity and acceleration 
of the mass, and the power requirement are low. There- 
fore, the use of the optimal control system in actual 
operations is feasible. Compared to the force control 
law method, detumbling can be achieved in one- fourth 
the time. This decrease is considerable since stabiliza- 
tion may require hours. It should be noted that the 
optimization technique need not only be considered from 
the standpoint of minimising time to detumble. Since 
time increases as mass decreases , a minimum mass solution 
can be obtained by fixing the time at the largest feasi- 
ble value. Viewing the comparison to the force control 
law method in relation to this examination of mass and 
time changes , it is possible that an object with a mass 
much smaller than that used in the force control law 
technique may be used to achieve simple spin in the same 
time period. When dealing with a mass of about 1,000 kg. 


-&JL, 




U & 


66 


this decrease will be quite considerable. This idea of 
mass decrease for equivalent times by using the optimiza- 
tion technique being investigated here can be applied to 
wobble damping discussed in Chapter IX. Since the 
angular velocities to be reduced are quite small, time 
may not be the critical variable. However, mass is 
always of importance in space applications due to cost 
per mass to be placed in orbit. The optimization tech- 
nique will reduce considerably the mass of the movable 
object needed to achieve simple spin in the same time 
period of active control. 

In regard to the nature of the optimal solutions that 
were obtained, the local minimum achieved here permits 
a faster detumbling time or a smaller mass when compared 
to existing solutions. In addition, the minimum seems 
to be an absolute minimum since the slopes of the and 
envelopes of oscillation for the 0.5% case shown in 
Figure 3 and the one percent case shown in Figure 7 are 
approximately constants. Specifically, noting that the 
constant slopes of these cases are comprised of about 100 
sec simulation time increments and that various guesses 
of the control variable were used, it becomes evident 
that there is a definite unique decrease in the peaks of 
the Wg an< ^ ^3 oscillations for each case. There seems 
to be only one minimum for a specific casej hence, it is 
an absolute minimum. 








67 




Considering actual operations in space, a hybrid 
computer may be more effective since modeling the vehicle 
dynamics on an analog will permit faster computation; as 
was discussed in the previous chapter, 100 sec of digital 
computer time are needed to obtain 100 sec of simulation 
control time. Comparison between computer predicted 
vehicle motions due to mass movement based on the optimiza- 
tion method and actual motions could then be continuously 
monitored, thereby permitting updating of the initial 
conditions for io^ 9 3 and for the subsequent simula- 
tion time increments for optimization. This updating 
could be done without having to stop active control. 

No attempt was made to improve on the first-order 
gradient solution; that is, an optimization method such 
as the neighboring extremal, which would utilize and 
require the solutions already obtained, was not used to 
further reduce the time or mass. If a digital system 
similar to the one used in this investigation was 
employed, the added computation time would necessitate 
periods of no active control since the total computer 
time would be greater than the simulation time. Even if 
a hybrid system could be used to sufficiently reduce the 
time needed for the added computation, the decrease in 
time or mass would probably be slight compared to the 
quantities already required. 
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The optimal control was applied specifically to 
stabilise a tumbling vehicle about its major principal 
axis. No computer runs were made to achieve simple spin 
about a specific geometric axis. However, the technique 
presented here can probably be applied to such a case by 
appropriately changing the performance index and the 
direction of mass motion- An indication of this is given 
by noting that the initial one percent mass case was 
actually made to spin about an axis approximately ten 
degrees from the maximum inertia axis of the whole system. 
This was due to the fact that the major principal axis of 
the main vehicle was chosen as the direction of the linear 
track and that the mass is large and far from the center 
of mass of the main body. In actual operations with 
large masses, the major principal axis of the system 
should be used as the direction for the linear track*, 
this was not done in this study since a close comparison 
to the smaller mass case was desired. To have spin about 
some geometric axes may, however, necessitate considerably 
higher energy input by the control system. This will 
probably occur for spin about axes close to the minimum 
inertia axis . 

It was assumed that the principal moments of inertia 
of the main body do net change. However, an explosion may 
result in part of the spacecraft blowing off. Also, 
tumbling itself may cause loss of part of the vehicle. 

- ■ A t-J 
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If this mass loss results in a large change in the 
moments of inertia, it should be included in the optimiza- 
tion. Specifically, the magnitudes of the moments of 
inertia should be corrected and the direction of the 
linear track should be altered to correspond with the new 
major principal axis if spin is desired about this axis. 

If direction change is not feasible, a redefinition of 
the performance index could be studied; but, as was 
shown, this would result in a longer time or a larger 
mass, and residual transverse angular velocity. 

Other types of control mass motions could be investi- 
gated using the optimization method presented here. 
Specifically, a translation other than linear or a 
rotation relative to the spacecraft . 

The optimal control technique investigated is signi- 
ficant in that it uses an open loop solution to control a 
vehicle in real time regardless of initial conditions. 

The highly nonlinea 1 equations of motion preclude the use 
of an optimal closed loop approach. Normally, a non 
optimal feedback method like the one proposed by Edwards 
would have to be used. The use of such optimal open loop 
solutions should be investigated for other active con- 
trol devices when optimization is desired. One area of 
application is for minimum time thrusting of control jets 
with constraints on three orthogonal components. 
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APPENDIX A 

EQUATIONS OF MOTION FOR A SPACECRAFT 
WITH A MOVABLE MASS 

The equations of motion are given below. The mass 
is permitted to move along the x axis which is fixed in 
the main body. 


w - [VAA(x 2 CMAl + x 4 CMA2 + CMA3) + 


VBBCx CMB1 + x 2 CMB2 + x 3 CMB3 + CMB4) + 


VCCCx CMC1 + x 2 CMC2 + x 3 CMC3 + CMC4) ] 


[1.0/ I A I ] 


= [VAACx CNA1 + x 2 CNA2 + x 3 CNA3 + CNA4 ) + 


VBB(x CNB1 + x CNB2 + CNB3 ) + 


VCCCx CNC1 + x 2 CNC2 + CNC3 ) H 1. 0/ | A| ] 


(I) = [VAACx CPA1 + x 2 CPA2 + x 3 CPA3 + CPA4) + 

z 


VBBCx CPB1 + x CPB2 + CPB3 ) + 


VCCCx CP^l + x 2 CPC2 + CPC3)Kl.O/ [A | 3 

where VAA = w 2 CWA1 + to 2 CWA2 + to w CCWA3 + x CWA4 ) + 
y z x y 


a) to CCWA5 + x CWA6) + w CWA7 + 
x z y z 


to x CWA8 + a) x CWA9 
y z 
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VBB = to 2 (CWB1 + x CWB2 ) + to 2 (CWB3 + x CWB4) + 
x z 

to to CWB5 + to to (CWBB + x 2 CWB7) + 
x y x z 

to to (CWB8 + x CWB9 ) + to xx CWB10 + 

y z- y 

x CWB11 

VCC = to 2 (CWC1 + x CWC2 ) + to 2 (CWC3 + x CWC4 ) + 

to to (CWC5 + x 2 CWC6 ) + to to CWC7 + 
x y x z 

to to (CWC8 + x CWC9 ) + to x£ CWC10 + 
y z z 

x CWC11 

|A| = x 4 u 2 CA3 + x 3 CA23 + x 2 uCA2 + *CA12 + CA1 
CA 3 = I x 

CA23 = - 2I ^ 2 y - 2i xzW 2 z 

CA2 = Vy + Vz - 4y - 4 + V** + V s * + 

i y vy 2 + i z pz 2 - 2i yz uyz 

CAX2 = -2I xy I ya yz - 2I xy u 2 yz 2 - 2I xy I z uy - 
2 I xyl i 2 y 3 - 2I yz I x2 uy - 2I xz y 2 y i z - 

2I xzV 7 ' " 2I xz li2z3 

- I I 2 - I 2 I - I 2 X - 

x yz xy z xz y 

» A. 



CA1 = I I I 
x y z 
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CWA1 

CWA2 

CWA3 

CWA4 

CWA5 

CWA6 

CWA7 

CWA8 

CWA9 

CWB1 

CWB2 


21 I i + I I yy 2 + I l„yz 2 - 
4 xy xz yz x x z 

2I x I yz wyz + I y I z ,jy2 + I y I z ,Jz2 + 


+ I y u 2 y 2 z 2 + I z u 2 y 2 z 2 + I z V 2 z 4 


^W 2 “ l 2 z - 2 “ 2 I yz ,l2y3z " 


2I yz ,j2yz3 ” I xy I xz wyZ “ I xy wy2 " 


I xz I xy’ jyz ' I xz wz2 


I + yyz 
y Z 


-X - yyz 
yz 


= I 


xz 


« yz 


= -i 


xy 


= -uy 


X y - I z - py 2 + vz 2 


= 2yy 


= 2yz 


= -I 


xz 


= -yz 








CWB3 = I 


xz 

CWB4 = yz 

CWB5 = -I - uyz 
yz 

CWB6 = I - I - yz 2 
z x 

CWB7 = y 
CWB8 = r 

xy 

CWB9 = yy 
CWB10 = -2y 
CWB11 = -yz 
CWCl = I xy 
CWC2 = yy 
CWC3 = -I 

xy 

CWC4 = -yy 

CWC5 = -I + I + yy 2 
y ^ 

CWC6 = -y 

CWC7 = I + yyz 
y 2 

xz 

. U & 


CWC8 = -I 


76 




CWC9 = -y 2 
CWC10 = -2 y 
CWC11 = + yy 

CMAl = l y y + I 7 y + y 2 y 2 + u 2 z 2 
CMA2 = y 2 

CMA3 = I y X 2 ♦ I y W 2 ♦ X a U2 2 - I 2 Z - 2I yz yyz 

cmbi = i yj7 yz * y 2 z 2 y + i z uy * u 2 y 3 

CMB2 = I u 
xy K 

CMB3 = y 2 y 

CMB * = Vxz + J xz^ z + I xy I z + I xy’ j y 2 

CMC1 = T yy * * -^y^ 2 * V 2 2 3 

CMC 2 = I y 
xz 

CMC 3 = y 2 z 

CMC 4 = I I + I yyz + I I + I yz 2 
xy yz xy HJ y xz xz^ 

CNA1 = I uz + y 2 yz 2 + I yy + y 2 y 3 
yz z 

CNA2 = I xy y 
CNA3 = y 2 y 






L A, . 


> 


C NA4 = I yz I xz + I xz py z ♦ I xy I z + I x y tJ y' 


CNB1 = 


-21 uz 
xz 


CNB2 = 


r . 2 2 

i x u + v- y 


CNB3 = i x i z + I x uy 2 + I z iiy 2 + u 2 y 4 + i z uz 2 + 

2 2 2 t 2 
^ y 2 “ X XZ 


cnci = i xz yy + J X y^ z 


CNC2 - 


2„„ 

u yz 


CNC3 = I xz I xy * Vyz + X *^ z * + 

y 2 y 3 z + I yz yz 2 + y 2 yz 3 

cpai = i y2 yy * vt 2 y 2 z + i y ys + y 2 z 3 


CPA2 = I xz u 


CPA3 = y z 


CPAU =11 + I uvz +11 + X uz 

urttH- X y yz xy My y xz xz H 


CPB1 = 


X xz^ + I *y WB 


CPB2 = 


y 2 yz 


CPB3 = I xz I xy * Vy, + ^ + I yz |jy2 + "V* + 
T 2,23 

i yz yz + u yz 


u A, ... 




CPC1 
CPC 2 
CPC 3 


-21 \iy 

2 2 

I x U + U z 

9 2 2 2 2 

i x i + i x yz + i y vy + u y ^ 


x 2 . ,2 4 t 2 

V z * W z " I xy- 
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APPENDIX B 
COMPUTER PROGRAM 


A listing of the first-order gradient optimization 
program is presented on the following pages. A flow 
chart is given in Figure 14. The following variables 
need to be specified for each case: 


EX = I 

x 

EY * I y 
EZ = I 

z 

EXY = I 

xy 

EXZ = I 

xz 

EYZ = 1 

yz 

X2MAX = “2 max 
X3MAX = “3 max 
X1INV = O3 x (t 0 ) 

X2INV = to (t Q ) 

X3INV = to z Ct 0 ) 

X4INV = x ( t q ) 

X5INV = &Ct Q ) = B<t 0 ) 






t. 
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YPO = y 
ZPO = z 
FORDIV = 

number of time steps 

IORDIV = 

TF = t f 
PEXLMH = x, 

PEXLML = x 1 
PECXVA = K 1 = K 2 

w = w 

EP = e 
SMAMAS = m 
BIGMAS = M 
YOU = y 
U = u x . 




•L & . 

















DIMENSION Ui 500) 

XENOY=*UOOE-20 

XENDZ=*1. G0E-20 

CX=, 6742093176*07 

£Y= ,62755678E*07 

E2=. 5152799 SE* 07 

EXY-0,0 

E XZ-O, 0 

E YZ-Oo 0 

EYX-EXY 

EZX-EXZ 

EZY-EYZ 

X2MAX=* 001 

X3MAX-o 001 

Sf3MC0A=0,0 

SOMC IB - 1,0 

$0MC2C = 0.0 

X5M A X- 1 <. 0 

XIiNV=. 103 

X2INV--,199 

X3INV=, 000286 

X4INV=OoO 

X5INV— 0,0 

P HNV-O, 0 

P2INV^0,0 

P3INV=0.0 

P4INV=0„0 

P SIN V— 0, 0 

WAN VA 1=10,0 

YPQ=5, 55 

ZPO=- 13, 7 

FORD! V=20o 0 

I0RD1V=20 

TF= 100.0 

PEXL IM= 2, 0 


MoosinisWd Moosuwirrv MoosiZNinv moositnifiv 4 toos w\ 

M005>£NJIIV MOOS>2NII3V MOOS JIN 1 1 IV M2iriV M2 4 : 
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0*GQX = V 
5*2-=! 
5 *2=H 


7PAF 2UI 5003 ? PAF3U* 500a, PAF5U1500) ? PAF6Ui500$, V*2J ? DElTAUJ5O0) 9 
8XST0R51 500) , PSTQR5! 5001 ? RSTOR5 (500 ?2 ) 

CA 1-E X*E Y^EZ-E X#E YZ^#2Hr Z^EXY^^-EY^EX 2 *0# EXY* EXZ* EYZ + 

IE X*E Y*YQU*YPO**2* E X*E ZnOU*Z PD**2-2 * 0*£X*EY Z*Y 0U* Y PO* ZPQ+EYfc EZ* 
2V0 : U*( YP0^2+ZP0^^2 ) *E Y^Y 01^*2 # ( Y PO* 4 +Y P0**2*Z°0** 2 ) +EZ* Y3U**2* 

31 YP3^2^ZP0^2->ZP0^4) ~EYZ*-2*YGU^lYPO««2<‘ZPO ; f £S 4 t 2 )-2»0*EYZ^Y0U** 2 
4*1 YPG**3*ZPG-EYPG*ZPG**3» -2 o 0*EXY*EXl *YOU*YPO*ZPO-EXY** 2^V0U^YP0 
5** 2- E XZ ** 2*YG U *2 P 0 **2 

CA2=E X*E Y-i-E X*EZ-£XY**2-EXZ**2 * 

lEX*YOU*YPO**2*EX*YOU*ZPa**2-»-EY*YCU*YPa**2*£Z*YOU*ZPO**2-2.G* EYZ 
2*Y0 U*YPO*ZPG 
CA 3=E X 

C A 2 3=~ 2 » 0 *E X Y* VO U * *2 * Y P 0-2 . 0 * E X l *' Y OU * *2 * 1 P C 

CA 12=-2.Q*EXY*{£YZ*YQU*ZPO+YOU**2*Y P0*2P0**2*EZ*Y GU* YPD*¥DU** 2 
1*YPG**3 )— 2* G*E XZ*! EYZ*Y0U#YP0*-Y0U**2*Y P0**2* ZPO^EY*YOlj*ZPQ+YOU 
2**2*ZPG**33 
CWA I=E YZ*YOU*YPO*Z PG 
CWA 2=— C WA 1 
CWA 3=EXZ 
CWA 4=Y0U*ZP0 
CWA5--EXY 
CWA6-- YGU^YPO 

C WA 7“E Y-E Z- YOU* YPG*YPO* YOU*Z PO*ZPO 

CWA 8=2o 0*YGU*YPO 

CWA 9=2* Q*YOU*ZPO 

CWB 1=-EXZ 

CWB 2~- YOU*ZPO 

CWB 3“E XZ 

CWB4=YGU*ZP0 

C WB 5--C WA I 

CWB6=EZ-EX-Y0U*ZP0*ZPG 
CWB 7= YOU 
CWB 8=E XY 
CWB9=Y0U*YPG 
CWB 10=- 2 o 0* YO U 


CWB11=-Y0U*ZP0 

CWC 1=£ XY 

CWC 2-YGUSYPO 

CWC 3=- E XY 

CWC4=-Y0U*YPG 

CWC 5=-E Y+E X* YOU*YPG*YPO 

CWC 6=- YOU 

C WC 7-C WA 1 

CWC 8=- E XZ 

CWC 9=- YOU-ZPO 

CWC XO-— 0#Y0U 

CWC 11=¥GU*YPG 

CMA1=E Y*YOU*FZ*YOU*YGU*YOU*lYPO*YPO*ZPO* ZPO! 

CMA2=Y0U*Y0U 

C MA 3=E Y^EZ+EY^YnU* YPG*YPO+ E Z*Y OU*ZPO* ZPO-EY Z* EY Z-2 .0* EYZ* YOU 
l^YPD^ZPO 

CMS 1=£ YZ*YOU*'ZPO* YOU*YDU*ZPG*Z PG*Y PO+-EZ*YOU*YPC+YOU* YOU* YPO** 3 
CMB 2-EXY#YGU 
CMB 3- YQU*YOU*YPG 

Cm 4-E YZ$E XZ*£ XZ*YGU*YPa*ZPO* EXY #£ Z*EXY *Y OU*Y PO*Y PG 

CMC 1=E YZ#YGU*YPG* YOU*YGU*YPO*YPO*ZP€H'EY *YGU* ZPO+YOU* YOU* ZPD** 3 

CMC 2=£XZ*YGU 

CMC3=YDU*YG0*ZPG 

CMC4-eXY*EY2*EXYSY0U*YP0*ZP0*EY*EXZ*EXZ*Y0U*ZPQ*ZP0 

CMA X=E YZ *¥Gti*ZPG* YOU*YGU*Y PO*Z PG*ZPO+ E Z* YOU*¥ PG*Y OU* YOU* YP 0** 3 

CM A 2~£ XY*YOU 

CNA 3=YGU*Y0U*¥PG 

CNA4-E YZ*E XZ*EXZ*YOU*YPO*Z PC.*- EXY*EZ*EXY *Y GU*Y PO*Y PG 

CNB1=~2-0*EXZ*YGU*ZPG 

CNB 2=EX*YQU*YGU*YGIJ-YPG*YPC 

CNB 3=E X*£Z*E X#YGU#YPO*YPO*E Z*Y GU*Y P0*Y PG+Y GU*Y OU*Y PG** 4 
1 *EZ * YG U*Z PO*ZP G* YG U*Y0U*YPO* Y PC*Z PO*Z PG- EX Z* EX Z 
CMC 1=E XZ*YQU*YPO*EXY*YGU*ZPD 
CMC 2=YOG*YDU*YPD*ZPG 

CMC 3-E XZ *E XY*E X*E YZ*E X*YQO*YPO*ZPG*EY Z*YGU*YPO*Y PO+¥OU*YOU* YPO 


1**3#ZPG+E yz*you*zpd*z po* YDU*Ya'J*Y PO*ZPQ**3 

CPA 1=E YZ *YG U^YPD-s- YGU*YGU*Y PC*Y PC*ZPO*EY *Y Of ZPG+Y OU* YOU* ZPO** 3 

CPA 2=E XZ*YGU 

CPA3=YOU*YOU*ZPO 

CPA4=EXY*EYZ*EXY*YOU*YPO*ZPO*EY*EXZ+EXZ*YOU*ZPO*ZPG 
CPB 1=CNC 1 
CPB 2=CNC 2 
CPB3=CNC3 

CPC 1=-* 2- 0*E XY*YQU* YPO 

CPC 2-E X*YGU* YOU*YOU*ZPG*ZPG 

CPC3=EX*EY*£X*YOU*ZPO*ZPG+EY*YOU*Y PO*YPO*Y GU*Y OU*YPG*YPO*ZPO*ZPO 
l-J*EY*YOU*ZPG#ZPG* YGU*YGU*ZP0**4-EXY*EXY 
24 T=O.G 

DEI TA=t IF- n /FORD I V 
DELTA2=DELTA/2 0 0 
WRfTEt 6,6101 

610 FORMAT? 8 IS 'TIME? SEC1 8 , 4 X 5 * WXl RAO/SEC 1 S2X , 5 WY I RAD/SECS * , 2X, 

l f, MZ CHAD /SEC I S2X, S MASSD3 SiFTJ S2X, 'DWX CRD/SC23 S 2X, 5 DWYIRD/SC2) fl 
2? 2X, ? DtfZ?RD/SC2J S2X, 9 DMASDSIF/$1 8 ? 2X , ®U CFT/SC2 1 9 3 
X? DaXUW-V 
X{ 21=X2INV 
XI 3 )-X3IMY 
XC4}=X4IMV 
X? 5 3=X5!NV 
X? 6 5=0. 0 
X? 7 1=0. 0 
XSTOR It 1 1=X? 1 ) 

XSTOR 2i 13=X? 2) 

XST0R3? 11 = X? 3) 

XSTOR 4? 11= XS 4} 

XSTOR 5? 13 =XC 53 
DO 2 K=X,IGRD1V 
DO 3 1=1,7 

R< I« U=F1U ,7»X,K,U,Y0U, PEXLIM j PECXV A, ci 

9PEXLMH, PEXCML , 



■ "X 





ft 



A 



ACAl,CA2»€A3fCA23,CA12 ?CVJA1 7 CWA2 7 CM A3 ,CW A4,CW A5, CW A 6 , CW A7, CW AS, 
BC WA 9 ? C WB 1 *C : fJB 2 *C WB3 ,C W84 * CWB5 , CVJB 6 *CW 97 * CW B 8 , CW B9 » CW B 10* CM B1 1, 
CCWCl,CMC2 s CWC3vCWC4vCWC5,CWC6 t CWC7,CWC8 *CWC 9 *CW CIO* CW Cl I? CM A In 
DCMA2 ? CWA3,CMBi ,CNB2 *CMB3 »CMB4« CMC1 *CMC2 , CMC3, CMC4, CNA1, CNA2? 
ECNA 3 9 CWA4 9 CNBl 9 CNB2?CNB3 9 CNCltCNC2< ) CNC3 9 CPAl 9 CPA2 9 CPA3,CPA4t 
FCPB1»CPB2»CPB3,CPC1 9 CPC2 ,CP€35 

3 XB 1 I) = Xt n+DELTA2*RU * 1 ) 

DXSTRXC K 5 =R(iai 
DX5TR2iK1=Rf 2*1) 

0 XS7R 3i K S “R 1 3* IS 
D XS7R 4t K ) =R ( 4 , 1 3 
TB= T+DELTA2 
T=T*DEtTA 
TIMEC 11=0,0 
TIME! K* 11 =T 
DO 4 1= 1*7 

Ri I*21=FH I f TBvXBjK 9 U,.V0U, PEXLIM»PECXV A, 

9P E XLMH ? P E XL M L , 

AC A 1 ?CA2 ?CA3 ,CA23,CA12 ?CWAL ,CMA2 9 CWA3,C^A4 P CWA5, CWA 6 , CW A7*CWA8, 
BCWA9,CWB1,CWB2 S CWB3»CWB4 vCWB5,CMB6 ,CWB7 , CWB8,CVJS9*CVi BIO, CWB11? 
CCMCl,eWC2»C1C3,CWC4 9 CWC5,CWC6 T CVlC7»CWC8 ,CWC9 ,€W CIO* CWC11, CM A 1* 
DCMA 2*CMA3 ,CMB1 *CMB2 ,CMB3 *CMB4 , CMC1 ,CMC2 ,CMC3, CMC4, CNA1, CN A2, 
ECNA3,CNA4 ? CNBl,CNB2tCNB3 9 CNCl ? CNC2 ,CNC3 f CP Al, CP A2, CP A3, CP A4, 
FCPB i 9 CPB2,CP83,CPCl,CPC2,CPC3) 

4 XB 1( I } = X( I 3*DELTA2#RtI ?2) 

DO 5 1 = 1,7 

Rt-I*3J=J?U£*TB,XB1 ,K T U»YOU, PEXL1M,PECXV A, 

9P E XLMH 7 PE XLML , 

ACA I-tCA 2,CA 3,C A23,C A12 ,CWA1 *CWA2 ,CWA3 , CW A4 »CW A5, CW A 6 , CW AT, CW A 8 , 
BCWA9*CMB1,CWB2,CWB3,CWB4,CWB5* CWB 6 7 CWB7 , CW 8 B , CW B9, CW810, CW 8 il, 
CCV}Cl,CMC2,CMC3,C WC4,CVJC5,CWC6,CWC7tCVJC8 ,CWC 9 ,CWC 10 , CWCU, CMA 1 , 
DCMA2,CMA3,CMBI*CMB2,CMB3yCMB4,CMCi,CMC2,CMC3,CMC4,CNAl,CNA2, 
ECNA3,CNA4,CNBl,CNB2*CNB3,CNCi ,CNC2 ,CNC3„ CPAl, CPA2, CP A3, CP A4, 
FCPBl,CPB2?CPB3vCPCl*CPC2,CPC3S 

5 XB( 1 )=XU S*DELTA*Rn ,35 


PEXLI M» PECXV A, 


DO 6 1=1,7 

6 R 1 1 * 4)=Flt I , T, XB , K ,U, YOU, 



9PEXLMH,PEXLML, 

ACAl ,CA2,C A3 ,CA23 ,CA12 ,CWA1 ,CWA2 ,CWA3 ,CW A4 , £W A5, CWA6, GW A7, CW AS, 
8C WA 9,C WB 1 ,C WB 2 » C W8 3 ,C W84 »CWB5 , C W86 ,C W B7 , CW B8 , CW B9 , CW B 10, CW 8 1 1* 
CCWCI,CWC2,CWC3 9 CWC4,CWC5,CWC6,CWC7,CWC8 ,CWC9,CWC10,CWC11,CMA1, 
DCMA 2,C MA 3 ,C MB 1 , CMB 2 ,C MB3 ,C MB4 « C MC 1 , CMC2 , CMC3, CMC4, CN A 1, CN A 2, 
ECNA 3,CNA4*CN8 1 ,CNB2 ,CNB3 ,CNC1 ,CNC2 P CNC3 , CPAl*CPA2* CP A3, CPA4, 
FCPB 1 ,C PB 2 ,C PB 3 ,C PC 1 ,C PC2 ,CPC33 
00 7 1=1,7 

7 Ml 5-XC 1 1’5'J DELTA/6® G3 R 1 1 ,1 J*2.0*CR( I ,2 3*RU,3 ) )*RU, 4) 3 
X5TORHK+13=Xm 
XST0R2C K*U =Xl 23 
XSTOR 3C K+li = Xt 33 
XST0R4C K*13=Xi43 
XSTOR5*K*l}=X{53 


2 CONTINUE 
PEN X71=X( 61 
PENXT2=XS 73 
K=XGRDI V+l 

TEMPORC 1 3 =XST0R1 1 1 0RDI V*1 3 
TEMPORf 23 = XSTGR2U OROI V*1 3 
TEMPORX 33 =XST0R3 U ORD I V+l 3 
TEMPORC 43 = XST0R4U8RDl V*13 
1=1 

DXSTR1* T0RDAD3 =F1(I ,T,TEMPOR,K,U,YOU, PEXL IM,PECXV A, 

9PEXLMH,PE XLML-, 

ACAI,CA2,CA3,CA23,CA12,CWA1 ,CWA2 , CWA3, CWA4,CW A5, CW A6, CHA7,CWA8, 
BC WA 9,C m 1 ,C WB 2 ,C WB 3 , C WB4 , CWB5 , CWB6 ,CW B7 , CW B8 , CW B9, OJ BIO, CW B 1 1, 
CCWC 1,CWC2tCWC3,CWC4,CWC5 v CWC6 ,CWC7 ,CWC8 ,CWC9,CWC1Qs CWCll* CM Al, 
DC M A 2 ,C HA 3 ,C MB 1 ,CMB2 ,C MB3 ,CMB4 ,CMC1 ,CMC2 ,CMC3,CMC4, CNA1, CNA2, 
ECNA 3,CNA 4,CN81 ,CNB2,CNB3 ,CNCl ,CNC2 ,CNC3 , CP Al , CPA2, CP A3, CP A4, 
FCPB1,CPB 2,CPB3,CPC 1 ,CPC2 ,C PC33 
1=2 


DXS7R2C I0R0AD3 =F1 (I ,T,TEMPOR ,K ,U ,Y CU, PEXL !M, PECXV A, 


9PEXLMH,PEXLML, 

ACAI,CA2,CA3,CA23,CA12,CWA1 ,CWA2 ,CWA3 »CW A4,CWA5 9 CWA6, CWA7,CWA8, 
BCWA9,CWB1,CWB2 V CWB3 S CWB4 ? CWB5»CWB6 9 CWB7 ,CWB8»CWB9, CWB1Q,CWB11 S 
CCWC1,CWC2 9 CWC3,CWC4,CWC5 S CWC6 9 CWC7,CWC8 ,CWC9 , CWCIO, CW Cl i, CMA 1* 
DCMA 2,CMA3,CMB 1?CMB2,CM83 ,CMB4 ? CMC1 ,CMC2 v CMC 3, CMC4, CNAUCNA2, 
ECNA3,CNA4*CNBi,CNB2,CNB3?CNCl ,CNC2 V CNC3 9 CPAU€PA2* CPA3*CPA4, 

FCPB 1 fCPB 2,C PB3 ,C PC 1 ,C PC2 9 CPC3) 

1=3 

DXSTR3C IORO AO 5 =FI 1 I 9 T 9 T£MPOR 9 K ¥ U 9 YOU, PEXl IM, PECXV A s 
9PEXLWH , PE XLML s 

ACA 1 ,C A 2 ,C A 3 ,C A23 *C AI2 9 C WAl 9 C WA2 9 CW A3 9 CW A4 9 CW A5 9 CW A6, CW A7* CW A8, 
BCWA9 9 CWBI 9 CWB2 9 CWB3 9 CWB4 9 CWB5 9 CWB6 9 CWB7 9 CWB8,,€W89 9 CWB10 9 CW Bll, 
CCWC 1,CWC 2?C WC3 9 C WC4 9 CWC5i>CWC6 9 CWCT ?€WC8 ,CWC9 9 CWClO, CKCllr CM Al, 
0CMA2,CMA3*CMBl 9 CMB2,CMB3vCMB4 9 CMClvCMC2 *CMC3 9 CMC4,CNA1 9 CNA2, 
ECNA3*CNA4*CNB1*CNB2,CNB3 » CNC1 9 CNC2 9 CNC3 , CP A1 , CPA2, CPA3, CPA4* 

FCPB If CPB 2 9 CPB3 *C PC 1 ,CPC2 ,C PC 3 ) 

1=4 

D XSTR4! I0RDA03 =FX 1 1 * T*T£MPOR ,K 9 U 9 Y0U 9 PEX t_ I M, P ECXV A» 
9PEXLMH,PEXLML 9 

ACA1vCA2,CA3,CA23,CA12 9 CWA1 ,CWA2,CWA3 ,CW A4,CW A5, CW A6, CW A7, CW A 8* 

BC WA 9,C WB 1 *C WB 2 9 C WB 3 V CWB4 ,CWB5 , CW86 , CW B7 , CWB8 , CW B9 , CW B 10, CW B 1 1, 
CCWC 1 *CWC 2<?C WC 3 9 CVJC4 9 C WC5 ?CWC6 , CWC7 ? CWC8 , CWC9 9 CWC10, CWCil? CHA 1, 
DCMA2,CMA3»CMB1,CM8 2 9 CMB3 ,CMB4 »CHCi ,CMC2,CMC3 9 CMC4 f CNAl,CNA2 9 
£CNA3tCNA4 9 CNBl 9 CNB2 9 CNB3 9 CNCl 9 CNC2 ,CNC3 1 CPAlt CPA2, CPA3, CPA4 V 
FCPB 1,CPB2,CPB3,CPC1,CPC2,CPC3 3 
DO 711 LK=1 9 TGRDAD 

WRITE(6,611) TI ME I LK) ? XSTORULKJ, XST0R2ILK), XST0R3CLK), 
1XST0R4U.K), DXSTRltLKI* DXSTR2 ILK) 9 0XSTR3CLK), XST0R5 ILK ) , UiLK) 

611 FORMAT! *0S10UXfE12.5}} 

711 CONTINUE 

WR I TE{ 6,611 3 TIME i IORDAD) ? XSTDRH IORDAD} 9 XSTQR2 UORDAD 3, 
1XST0R3! IORDAD) , XSTGR4TI DRDAD) , 0XSTR1 UORDAD) 9 DXST R2t IORDAD) , 
2DXSTR3I IORDAD) , XST0R5 U ORD AD) , UUGR0AD9 
WRITE (6,612) TF 

612 FORMAT! »0»,F15.6) 


WRITES 6*7005) PEN! T1 , PENI T2 
7005 FOR MATt , 0 , t2{3X»Ell«4|}} 

00 4441 KIK=1,I GRD AO 

FORCEX-YOU# ! UC KI K) +YP0*XSTDR1 { KIK) *XSTDR2lKl K3-XSTOR4 IKXK )* 
lXST0R2tKIK)*XSTaR2iKI K3 -XST0R4 SKI K3 #XST0R3 {KIK3*X$TOR3lK IK J 
2*ZPO*XSTOR1{KIK) #XS70R3 ( KI K) P0^0XSTR2 l KIKS— YP0*0XSTR3IK IK 5 ) 

TDO T“FGRCE X*XSTGR5 { KI K3 

R IN XDN=XSTOR 5C KI K J *ZPG*XST0R2 ( KI K3 -Y PG*XST 0R3 t K IK 3 
R IN YDN=XSTGR4l KI KJ *X5TGR3 {KX K3 -ZP0*XST0R1 l KIK3 
R INZDN=YP0*XST0R1 { KI K) -XSTQR4 { KI KJ *XST0R2 ( KI KI 

ENEKIN=* 5#{ E X^XSTORl { KI K) *XST0R1 IKIK) +EY«XST0R2 IK IK )*XST0R2CK IK ) + 
1EZ#XST0R3{KIK)*XSTGR3{ KI K) )+.5*Y0U*{RINXDN*RINXDN*RINYDN#R INYDN + 
2R IN2DN*R INZDN1 -XSTORl t KI K) *X ST0R3 1 KI K} *£XZ-XST 0R2 C K IK )* 

3XST0R3I KIK)*EYZ— X5T0R1IKIK1 *XST0R2 { KI K3 *EXY 
ANGMO X~E X*XST0R1 t KI K) -£ XY*XST0R2 t KI KI -EX Z*XST 0R3I K IK } +YOU* { YPG* 

1R INZDN-ZPQ^RI NYDN3 

ANGMOY=-EXY*XSTORi C KI K) *EY£X$T0R2 ( KI KI -£Y Z*XST OR 3 { K IK 3 *YQU* { ZPO* 
IRINXDN- XST0R4I KI K3 *RI NZDN) 

ANGMOZ-— E XZ *XSTQRl { KI K) -EYZ*XST0R2 1 KI K 3 *E Z*XST0R3 < K IK 3 *YGU* { 
1X5T0R 4< K1K?*RINY0N-YPQ*RINXDN) 

ANGMOT= SORTS ANGHQ X*ANGMOX*ANGMOY*ANGMOY * ANGMOZ* ANGMOZ 3 
WRITE! 6,4442! TIME t KI K3 ,FGRCEX ,TDOT vENEKIN, ANGMGT, ANGMGX, ANGMQY, 
1ANGMQZ 

4442 FOR MATt , G , ,8SlX,Ei2,5) ) 

4441 CONTINUE 
TP^TF 
TFP=T 

PDELTA=-DELTA 
P0ELT2=P0ELTA/2.0 
DO 14 M=l,3 
GO TO C 15,16,20) , M 
15 PUt-PlINV 
PS 2 3=P 21 NV 
P{3}=P3INV 
PS 4)~P4INV 


PI 5?=P5INV 

P STOR 1 U GRD I V* 1) Cl 3 
P STOR 21 IORD I V* 1 3 =P C 2 1 
PSTQR31 IORDI V+l ) -PC3^ 
PST0R4S IORDI V*I> -Pi 4? 

P STOR 51 IORDI V*XJ =P( 5} 
SOMTCG= SOMCOA 
SOMTC 1- SOHC IB 
SOMTC 2- SOMC 2C 
GO TO 17 
16 PI D^OeO 
Pi 2 !=Ci,€ 

Pi 3 )=0 0 0 
P(4I-0*0 
PC 5)=0<,0 

RSTGRH IORDIV*L*l) =PUJ 
RSTOR 2i IORDI V+I,U =Pt2? 
RST0R3C I0RDXV*-1»U=PI3> 
RST0R4C IORDI V-J-lsl) ^=PI4J 
RSTOR 51 IORDI V*l,l) =PC5) 
SOHTCQ=0*0 
SOMTC I=Qi 0 
SOMTC 2=0<» 0 
GO TO 17 
20 PC 1 ^0 o 0 
PC 2H.0 
Pi 35=Oi 0 
P{ 4>-0i 0 
Pi 55=0, 0 

R STOR IS IORD IV* 1*2) =Pill 
RST0R2C IORDI V*1 ?25 —PS 21 
RST0R3I IGRDIV*1,2? =Pi3) 
RSTOR 41 IORD IV* 1,2? =P*4? 
RSTOR 5( IORDI V*1 ,2) =PI5) 
SOMTCO^Oo 0 


SBMTC 1=0 oO 
SQM7C 2=0o 0 
17 DO 8 J=l* SORDI V 
OP 9 1=1,5 

RPC Lt 1} =F2I L,TP,P, SOHTCOs Jt U 9 yOU»XSTOR5i M, 

9PEXLMH 9 P£ XL ML 9 

AC A 1,CA 2 ,CA 3 sCA 23 ,C A12 ? CWA1 , CWA2 »CWA3 9 CWA4, CW A5* CWA6? CWA7* CWA8, 
BCWA9,CWB1 9 CWB2 9 CWB3,CW84,CWB5 9 CWB6,CVIB7 9 CWB8 9 CW89 9 CWB10, CWB11, 
CC WC 1 ,C WC 2 ,C WC 3 ,C WC 4 ,C WC5 ,C WC6 , C WC7 , CW C8 7 CW C9 , CW CIO, CW C 11, CM A 1, 
OCMA 2,CMA3,CMB 1 ,CMB2 ,CMB3 9 CMB4,CMC1 9 CMC2 ,CMC3,CMC4 9 CNA1, CNA2 9 
ECWA3 jCMA4,CNBl ,CNB2 ,CNB3 ,CNCL ,CNC2 ,CNC3, CPA1, CPA2, C^A3,CPA4, 
FGPBi,CPB2,CPB3,CPCl,CPC2 9 CPC3 9 
1SGMTC1, 

2SQMTC2, 

4XST0R1, XST0R2, XST0R3 9 XST0R4, IQRDIV, PEXLIM, PECXVAJ 
9 PBi L 5-P C L)+PDELT2^RPC 1*13 
TPB=TP*PDELT2 
TP=7P*PDELTA 
OQ 10 1=1,5 

RPf L r 2?=F2n_,IP8,PB,S0MTC0,J,U,Y0U,XSTQR5T M, 

9PEXLMH,PEXtML, 

ACA1 *CA 2, C A 3, C A 23, C A 12 ,CWA1 ,CWA2,CWA3 *CW A4 ,CW A3* CW A6 9 CW A7, CW A8 S 
BC WA 9 ,C WB 1 ,C WB 2 , C WB3 , C WB4 , CWS5 , CWB6 ,CW B7 * CW B8 , CW B9 , CW BIG, CW B 1 1, 
CCWC 1,CWC2,C WC3,CWC4,CWC5 ,CWC6 ,CWC7 *CWC8 ,CWC9 , CW CIO, CW Cl I, CM A 1, 
OCMA 2,CMA3,C MB1 ,CMB2 ,CMB3 9 CMB4,CMC1 *CMC2,CMC3 9 CMC4, CNAI?CNA2, 
ECNA3,CNA4,CNB l<7CNB2,CNB3,CNCl,C!viC2 ,CNC3 9 CPAl, CPA2, CP A3? CPA4 9 
FCPB 1,CPB2,CPB 3,CPC 1 9 CPC2 ,C PC3 , 

1S0MTC 1? 

2S0MTC2, 

4XST0R1, XST0R2, XST0R3 , XST0R4, IOROIV 9 PEXLIM, PECXVAJ 
10 PBl€U=PiU+PDELT2*RPiL,2} 

DO 11 1=1,3 

RPH:?33=F2{1_,TPB,PB1,$OMTC0 7 J,U,YOU ? XST0R5 , M s 

9PEXLMH*PEXLML, 

ACAl,CA2,CA3,CA23*CAl2,CWAi,CWA2,CWA3*CWA4,CWA5,CWA6* CWA7,CWA8, 


BCWA9*CWB1,CWB2,CWB3 »CWB4sCWB5 ,CWB6 ,CWB7 * CWB8 y CV? B9, CW 810, CVsfBll, 
CC WC 1 *C WC 2 ,C WC 3 *C WC 4 f C WC5 v C WC6 7 C WC7 9 CWC8 , CW C9 7 CW CIO , CW Cl 1, CM A 1, 
DCMA2*CMA3,CMB1 sCMB 2 ,CMB3 sC MB4 sCMCl *CMC2 ,CMC3,CMC4,CNA U CNA2 5 
ECNA3?CNA4,CNBl T CNB2,CNB3sCNCl , CNC2 , CNC3 s CPA1 , CPA2, CPA3, CPA4* 
FCPB 1 sCP 82 9 CP83 *CPC 1 ?C PC2 jCPC3 7 
ISOM TCI, 

2SDMTC2, 

4XST0R1, XST0R2 * XST0R3 , XST0R4 9 TORDIV* PEXLIM, PECXVA) 

11 PB{ L J=P{ L)-*-PDE LTA^RPC L,33 
DO 12 L~ls 5 

12 RPi Lt4)=F2{LsTP,PB ,S0MTC0*JsU,Y0UsXSTaR5 „ M* 

9P E XLMH y PE XLML 9 

ACAlsCA2sCA3*CA23sCA12fCWAl ,C UA2 ,CWA3 T CW A4 » CW A5* CW A6, CW A7,CWA8, 
BC WA 9, C WB 1 ,C WB 2 * C WB3 * CWB4 * CWB5 , CWB6 9 CW B7 , CW B8 , CW B9 , CW B 10, CW B1 U 
CC WC 1 9C WC 2cC WC 3 ,C WC4,C WC5 S CWC6 sCWC7 sCWCB t CWC9 , CWC10, CWCU9 CM AU 
DCMA 2 ?CMA3 »C MB1 9C MB2 9 C MB3 «CMB4 , CMC1 t CMC2 *CMC3s CMC4, CNA1, CNA2, 
ECNA 3,CMA4jrCN8 l,CNB2vCWB3 ?CNC1 ,CNC2 P CNC3 , CP A1 , CP A2? CP A3, CPA4« 
FCPB lsCPB2?CPB3 ,CPC 1 sCPC2 9 CPC3 9 
1S0MTC 1* 

2S0MTC 2s 

4XS7GRI* XSTQR2v X5TQR3 * XST0R4, IGRDIV* PEXLIM* PECXVA) 

DO 13 L=l*5 

13 Pit 1-P< U-HPDELTA/BcO) *{ RP{ L* U*2 oO*{ RP ( L* 2 3 *RP (Ls 3 H +RP *L , 4) ) 
GO TO I 18s 19 ,21 ) sM 

18 P STORK IORDI V*l— J) =P( 1) 

PSTGR2I IORDI V*l-J3 =P?2> 

P$T0R3< IORDI V+l-J) =P 13 i 
PST0R4I IORDI V+l- J) =P*4) 

PST0R5I IORDI V* 1— J) =P(5) 

GO TO 8 

19 R STORK IORDIV^l-0 ? l)=Pil) 

R STOR 2( IORDIV^-1-Jsl) =PI2) 

RST0R3C IORDI V+l-J , 1 ) =P { 3) 

R STOR 41 IORDI V*l— J 9 1 V ~Pt43 
RST0R5I IORDI V*l-J,15 -P(5) 


r 


GO TO 8 

21 R STORK IORDI Vr 1— J*2) =Pf 19 
RSTOR2! 10RD I V-M—J s2S =P(2S 
R STOR 3( 10RDIV-M-J*23=PC31 
RST0R4C IORDIV*!- J,2§ =Pf4} 

R STOR 51 SGRDIV*l~J*2)=P!5i 
8 CONTINUE 
14 CONTINUE 

DO 22 K-UIORDAD 

C0EFAS=eAI*CA24XSTGR4f K? **2*Y0U*CA3*XST0R4(K)**4*Y0U**= 2 
l-f-XST0R4CK)4CA12*XST0R4(K) 4*3 4CA23 
IF C A8 5! Ul KS ). LTcPEULI Ml GO TO 33 
IF (UfKI.GT.0c03 GO TO 34 
l IF (Uf KJ.LTc0.03 GO TO 35 

K 33 PeNC02=0o0 
y GO TO 36 

34 PENCG2=PECUVA 
GO TO 36 

r 35 PENCG 2=~PEC OVA 

36 PAF 1U(K> =! 1. O/COEF AS) 4(CViBl 14 ( XST0R4 (KS 4CMB1 +XST0R4( K >4XST0R4(K) 

^ 1*CHB 2+XSTDR4! K5 4*34CMB3*CMB4 ) *CViCll*iXSTGR4 ! KS* CMC1*XST0R4 
2( K 3 4 XSTOR 4! K ) 4C MC2 ^ X ST0R4 ( K J 4^3 *C MC3-S'CMC4 } ) 

PAF 2U( K S =i 1« O/COEF ASS 4(CVJB1 1 4(XSTGR4 ! K5 4CN81+XST 0R4CK ?4XST0R4{K } 

I 14CN8 2*CNB33 *CWC114(XST0R4f KJ 4CrCl*XST0R4(KS4XSTGR4{K34CNC2*CNC3t J 

\ PAF3UC K S~C 1. O/COEF AS1 4( CVJB1 1* C XST0R4 (KS 4CP81-5-XSTOR4C K S4XSTGR4(K )4 

2CPB 2*EP83I+CWC U4f XST0R4 ( KS 4CPC1+XST0R4 « Ki 4XSTQR4SK 84CPC2*CPC3S } 
PAF5UC K 0 

! PAF6U(K5-2.04PENC024(ABS(U(KJ S-PEULJM3 

? AIHNIEKS-URSTGRI fK t i)4PAFlU(K)*RSTQR2(K ? 1)4PAF2UUK)+RST0R3{K 9 1) 

' 14PAF 3UC K )*R STOR5I K *1 J 4PAF5 U f K) * PAF6U { K) ) 4#2 )/W 

i A I f IN2f K) =1 RSTQR1 ( K*1 ) 4PAF1 Uf KS + RSTOR2 <K T lS4PAF2lHKS frRST DR3fK* 1) 

14PA F3U( K S *R ST0R5C K ?1 S 4PAF5 U( K) *PAF6U(K) ) 4 (RSI ORH K, 2 1* PAF 1U (K 
2r+RST0R2«K,2S*PAF2U( K)+RSTOR3(Kv2J4PAF3U nO+RSTOR5(K, 2J4PAF5UCKS 
3*PAF 6Ut K S J /N 

I AH 1N4( K) - t ( RSTOR1 ( K »2 J 4PAF1 U! KS +RSTOR2 ( 2J4PAF2U € K ) *RST0R3(K* 2)4 

k 



I If* — 


« 1PAF3U! K)*R$7GR5{ K*2) #PAF5U{ K) + PAF6U1K) J**23/W 

AIJINim=(PSTQRll KWAFlUUUFPSTOR2m*PAF2UIK3*PSTOR3&K3* 

2PAF3UI KJ-i-PST0R5f K3 *PAF5U( K) *PAF6Ui K} I *C RSTGRHK? 1 )* PAFIU IK J 
3+R STOR 21 K, 1 3 *PAF2 U C K3 * RStOR3 l K ? 1 3 *PAF3U I K J *RS T 0R5 1 K? 1 3*P AF 5U * K ) 
4*PAF6UCKH/« 

AIJIW2I K3=C PSTORHK3 *PAF1 U? K) + PSTOR2 ( K3 *PAF2U IK 3 S-PSTGR3IK 3*P AF3U! < 
1 3 -H* STOR 5$ K ? =fcPAF5UI KHPAF6UC K) } *1 RSTORU K*2 »*PAF1U IK3 *RSTGR2IK? 2 
2 W AF2U3 K 3 *RST0R3? K?23 *PAF3UI K3 *RSTOR5 t K ?2 J*PAF5U C K J *PAF6U IK 3 ) /H 
AJ^INTC 5 PSTOR1 I K) *PAFlUi K3 *PSTOR2 (K3#PAF2UlfO+PSTOR3CK J*PAF3U 
2iK?*PSTGR5S KJ*PAF5UIK3*PAF6UI K3 3 #*2J/W 
22 CONTINUE 

AIT! I? 1 ?=FNTGRL { I GRDi V+l ? DELTA? AIIIN1J 
A Ilf 1„23=FN7GRL! IORDI V*1 , DELTA? AIIIN23 
AIH2,13=Aim,23 

AIfC2?23=FNTGRLMGRDIV*l, DELTA? AIIIN4 3 
V AIJCH =FNTGRLMORDI V*l? DELTA? AIJINXI 

K AUf23 =FNTGRLf IGRDIV+X , DELTA? AIJIN23 

y AJJ =FNTGRUI0RDIV*1, DELTA? AJJINT3 

HRITE(6?6143 AIHL«U« AXIU*23? AII{2 S 23? A1J113? AiJ 12 3, AJJ 
614 FORMAT! ®G* 5X?£1 Q 0 3U 
r* BETMV=Aim,X3*AI *12 *2) -AH 11,23**2 

DLPSX1=-EP*PENIT1 
^ DLP SI2-— EP^PENI T2 

Vf 1 i=~i l^G/BETTO*!! All C2, 231 *f DIPS Il*AIJU3 U(-AII(l»2i 3* CDLPSI2 

V€ 2 1.0/DETtiV}*U-*AXH2 v l} }*{DLPSI 1*AI JUJ 3*i AIf{ 1? 13 3* CDLPSI2* 

1 iAij( 2 m 

IFf A IH 1? 1? «GTo 0* 03 GO TO 7701 
GO TO 7703 

7701 IFI A I If 2, 23 «GT« 0. 0) GO TO 7702 
r DETM V=A I H 1 ? 1 3 

DLP SI I=-EP*PENX T1 

r ot p si2-0e0 

[ V* n=-C 1® O/DETMV? *IDLPSI1+ATJC133 

[ Vt 23=0.0 


(D 

C71 



GO TO 7702 

7703 IFCAIH 2*2*. GT„ 0*01 GO TO 7704 
DETMV=0*0 

DLP 51 1-0*0 
DLP SI 2=0*0 
VC 1NG.0 
VC 2 NO* 0 
GO TO 7702 

7704 DE7MV=A 11(2*23 
DLP SI 1-0*0 
DLPSI2=-EP*PENIT2 
VC 13=0.0 

VC 2)=-C X.O/DETHV* #CDLPSI2*AI JC23 3 
7702 CONTINUE 

DETI S=AI Id *13 *AM (2 *21 -All 11*23**2 

GOTVAI=AJJ-f 1-0/DE TIS>*UAIJU I *Ai 1 12,21 -AIJC2NAIH 2*1 n*AIjm + 
11 -A um*A-ZItl?2!*AI JC25*AII [1 *151 *AI JC2 33 
WRITE! 6 5 361 3 DETMV* Vd) » VC23 
361 FORMAT! *0* *4C3X*E10®333 

IFtABSC XSTQR2C IORDIV-M.il . LT*XENDY3 GO TO 25 
GO TO 29 

25 IFCABSC XST0R3C XORDIV^in.LT. XENDZ? GO TO 27 
GO TO 29 

27 IFC ABSCSOTVAH .LT. WANVAU GO TO 28 
29 CONTINUE 

DO 23 K=1 T I BROAD 

DEL TAUC K NPAFH3I K3 H PST0R1 1 K) + RST0R1 IK* 1 3*V i 1 3 *KST0fUC K* 23*V 1 2 1 3 
l*pAF2Uf K**l PST0R2C K3*R5TGR2 CK*1**VUJ +RSTGR2CK* 2 )*V ! 2 3 3*PAF3UK) 
2»5P STOR3C K1*RST0R3(K*1!*V< 1H-RSTQR3 ifc*2 3 *V (2) NPAF5U CK 3*{P$TGR5?K 3 
33-RSTOR 51 K* 1 3 *VC 1 J +RST0R5 C K*2 1 *V 1 1 1 ) +PAF&U ( K) * IV U ) +V < 2 H 
Ui K }=UC K3— DELTA UC K3 /W 
23 CONTINUE 
GO TO 24 

28 RETURN 
END 


FUNCTION Pin *T*X 9 K,U»Y 0 U, PEXLIM,PECXVA, 

9P E XLMH , PE XL ML p 

ACAI^CA2^CA3 ff CA23vCAi2»CWAl , CHA2 ,CUA3 , CW A4,CW A5, CWA6, CWA7, CW A8, 

BC VJA 9 ,C m l ,C m 2 ,C WB3 , C WB4 , CWB5 , C WB 6 , CW B7 , CW B 8 , CW B9, CW BIO, CW B i 1, 
CCWCl,CWC2,CVJC3,CiC4,CV3C5,CWC& ¥ CWC7 s CWC8 9 CWC9 9 CWCI0, CWCIUCMAl* 
DCMA2,CMA3,CMB1,CMB2 9 CMB3 ,CM94 ,CMC1 ,CMC2 , CMC3 , CMC4, CNA 1 , CNA2, 
ECNA3,CNA4 9 CNB1 9CNB2,CNB3 , CNCI ,CNC2 ,CNC3 9 CPA1 , CP A2, CP A3, CPA4, 

FCPB IpCPB 2,C PB3,CPC 1 ?CPC2 ,CPC31 
DIMENSION X£ 7) , UC 5001 

COEF AS=CA1*CA2*X*4J *#2#YOU*C A3*X 14} **4*Y0U**2*X f 4J*C A 12*XI4)** 3* 
1CA23 

VAA^XI 21*XS 2)*CVJAl*XI3>*Xt31 *CWA2*XU i*X «2mCW A3+XI41*CWA41 
1-i-XI l)*XS3WCVJA5*X(4] *CWA61 *X { 21^X13) *CWA7*X {2)*X (5)*CV1A8*X( 3} 
2*X{ 5 1 *C WA 9 

VBB"XI 1>*XC l)*tCWBl+XC4) *CWB2 1 *X 13 ) *X m*£CWB3*X{41*CWB41 
1*XI 1)*XI 21*CWB5**X£ 1J *Xm *(CWB6*X«4)*X?4 &*CWB71+Xt 2S*XC3$* ICWB8 
2+XI 4)*CWB91 + XI21*{ X54J *X£5$ *CWB10 3 I K3 *CW Bl i 
VCC = Xm*xm**CWCl*Xt4?*CWC2)*X{2>*X{2 5*ICWC34'X{4l*CWC4HX( 1 5* 
lX<2mCNC5*XI4S*Xl4)*CWC6)*XU}*Xl35*CWC7*X<2l*xm*iCWt8*X(4)* 
2CWC9)*X( 3)*« X145*X«5> *CWC10mJCK3 *CWCU 
VC M A = XI 4 $ * X 1 4 } M A l* X 4 4 1 **4 *C M A2 * C M A3 
VCMB— XI 45SCMB1* XI 4}*XI41 #CMP24-X?4I **3 *CMB3+CNB4 
VCMC = XI 4mMCl*XI41*XI41*Cno2-a-X{45**3*CMC3'5-CMC4 
VCNA-XC 41 $€NAI*Xt 45*XI41 #CNA2*X 141 **3 *CNA3+CNA4 
VCNS-Xf 4) *.C-NB1*XI 41 *XC41 *CNB2+CNB3 
VCNC-Xf 4) #CNCI*: XI 45 $X441 4CNC2-S-CNC3 
VCPA = XI 4}*CPATS:XI4i£XI41*CPA2*Xi41**3 *GPA3+CPA4 
VCPB-XI 4) *CPS 1+XI 41 *XI4$ *CPB2*CPB3 
VC?C=XC 4) PC 1*X14)*XI4$ *CPC2*CPC3 
GO TO C 1 s 293»495 9 6, 7S 9 I 

1 Fi=4 IcO /COEFA SI *£ VAA*VCMA* VBB*VCMB4-VCC*VCN]C) 

RETURN 

2 F 1=1 UO/COEFASI*! VAA^VCNA* VB8*VCNB+VCC*V CNC5 
RETURN 

3 F1-! UO/COEFAS)*! VAA*VCPA+ VBB*VCPB+VCC*VCPU 


r 








RETURN 

4 F XI 5 3 
RETURN 

5 F 1“ Ui K 3 
RETURN 

6 IFI* X{ 41 «LTePEXLMHI GO TO 7001 
F I-PEC XVA*{ X( 4) -PE XLMHsl 

GO TO 7002 

7001 F 1- 0 » 0 

7002 CONTINUE 
RETURN 

7 IFi X4 4}»GT„PEXLMU GO TO 7003 
Fl=PECXVA*l-X{ 45-5- PE XL MU 

, GO TO 7004 

\ 7003 F 1= OoO 
i) 7004 CONTINUE 
RETURN 
END 

r 

FUNCTION F2I LvTP v P ? SONTCOv JsUvYOUsXSTORS^ M* 

9PEXLMH-pPE XL ML 9 

ACA 1?CA 2«CA3?CA23 ?CA12 9 CWAI 9 CWA2 ,CWA3 ,CW A4 9 CWA5 9 Ctf A 6 * CW A7 9 CUA8 9 
BCUA9? CWB 1 9 C NB 2 ?C WB3 9 CWB4 9 CNB5 ? CWB 6 e CWB7 9 CN 889 CWB 9 V CWBIO 9 CWBII 9 
CC WC 1 9 C WC 2 v C WC 3 ,C UC 4 „C V3C5 9 C WC6 9 C WC7 9 CW C8 9 CVJ C9 , CW C 10 9 C W C 1 19 CM A U 
C DCHA29CNA39CMBI9CNB29CMB3 9CMB49CMCI 9CMC29CMC39CMC49CNAI9CNA29 

2 CNA 39 CNA 49 CNBI 9 CNB 29 C NB3 9 C NCI 9 CNC 2 ?CNC3 9 CPAIv CP A2v CP A3* CPA4 ? 
j FCPB I 9 CPB 29 CPB 39 CPCI 9 CPG 2 ?CPC3 , 

r ISONTCI 9 

2S0M7C29 

t 4 XSTORI 9 XST 0 R 29 XST0R3 9 XST0R4, IORDIV 9 PEXLIM 9 PECXVA) 

i DIMENSION P*5?, U« 500 ) 9 XSTORU500S 9 XSTaR2(500| v XSTOR3<500)9 

1XST0R4C 500} 9 XSTGR5C500) 

JJ“IGRDIV*2-J 

| VCMA=XST0R4S JJ3*XSTGR4l JJS *CMA1*XST0R4I JJ) **4*CMA2+CMA3 

k 



i 


-is a 


VCMB-XST0R41 J J ) MB1* XSTGR4 ( JJ J *X$TGR4 U J? *CMB2*XST0R4t J J I** 3 
1*CMB3*C MB 4 

VCMOXSTOR4* J J> *C MC1+XSTGR4 1 J J? #XST0R4{ J JJ #CMC2 +XSTGR4C J J 3 

l*Ch*C3*CMC4 

VCNA=XST0R4{ JJ) *CNA1*XSTGR4 ( J J ? *XST0R4 1 J J) *CNA2*XST0R4 U J 3** 3 
1#CNA3*CNA4 

VCNB=XST0R4! JJ? *CNBl*XSTGR4 U J) *XSTGR4 [ J J ) *CNB2*CNB3 
VCNC-X5T0R4* JJ? *CNC1* XST0R4 U J ) *XSTGR4 (JJ) *CNC2*CNC3 
VCP A=XST0R4{ J J? ’CPAI+XST0R4C JJ?*XST0R4< JJ) *CPA2*X$T0R4(J J )**3 

1#CPA3*CPA4 

VCPB=XSTGR4{ JJ)*CPBi*XST0R4( JJ?*XST0R4( JJ?*CP82*CPB3 
VCPC-XSTGR4C J Jl ^ PCI* XSTGR4 i J J? *XSTGR4 ( J J) *CPC2*CPC3 
CGEFAS=CAi*CA2*XSTGR4{ JJ? **2 *YGU+CA3*XSTGR4{ J J3 **4*Y0U** 2 
1*X$TGR4{ J J ? #CA 12* XST0R4 ! J J) **3 *C A23 
GO TO l 1s 2 9 3*4*5) s L 

1 P AF 1X1" l UO/COEFASS*! VC MA X ST GR2 1 J J? *« CW A3 *XSTGR4{ J J ?*CW A4J 
l*XSTGR3t JJ?*iCWA5*XSTGR4UJ?*CWA&n*VCMB*£2oG*XSTGRlSJjmcWBl 
2-5-XSTGR4I J J? *CWB2) * XSTQR2 ! J J? *CWB5*XST 0R3 I J J! * * CW B6*XST0R4£ J J 5* 

■ 3XST0R41 JJ)*CWB7?I*VCMC*{2«0*XST0RH J J ?* C CWC1*XST0R4S J J ?*CWC2> 
4*XS TDR 21 J J ? *t C VC 5* XST0R4 ( J J? *X STQR4 ( J J5 *CW C6 ) *XST GR3 {JJ )£CWC7) ) 

, PAF2Xl=£loG/CGEFAS3*( VCNA*£ XST0R2UJ) * E CW A3*XST GR4U J )*CWA4) 

1 *XSTGR 3{ J J ) *{ C WA S* XS T0R4 C J J 3 *C WA6 } ) +V CN B* 1 2 o 0*XS T OR 1 ! J J ? * i CW B1 
2*XSTGR4{ JJ? *CWB25* XST0R2 l JJ? *CWB5*X$T0R3 C JJ?*lCWB6*XST0R4i JJ )* 
3XSTQR4S JJ \ *C WB7S ? * VC NC *1 2* 0*X S TORI ( J J l * £ CW Cl *XST 0R4 { JJ 1*CWG2! 
4*XSTGR2I JJ? *£CWC5*XST0R4{ JJ) *XST0R4 ( J J3 *CWC6 )*XSTOR3( J J 3*CWC73 ) 
PAF3X1-C lo O/COEF AS) VCPA^t XSTGR2 { JJ? *{CW A3*XST 0R4CJ J ?*CWA4J 
l*XSTOR3{ JJJ*tCWA5*XST0R41 JJ? *CWA6n*VCPB*(2.>0*XST0RHJJ ?MCVJB1 
2+XSTGR4? J J) *CWB23 * XSTGR2 { JJ) *CWB5*XSTOR3 ( J J3 *t CWB6*XST0R4{ J J 3* 
3XSTGR4i JJ3*CW8 73 )* VCPC*{2« 0*XSTGR1 l JJ )* { CWC1 *XSTGR44 J J )*CWC2> 
4*XSTGR2{ JJ)*{CWC5*XST0R4{ JJ3*X$TOR4( J J? *CW C6 3+XSTGR3* J J 3*CWC7) i 
F 2-— P i 13#PAFlXl-Pl2>*P4F2X!-Pt3>*PAF3Xl 
RETURN 

2 PAF IX 2= { lo O/CQEF A S3 *{ VCMA*E2 o 0*XST0R2 i JJ)*CWAl*XSTORHJJ )* { 

1C WA3*X STORES J J) *CWA4) *XSTOR3 { JJ? *CWA7 *XSTGR5 { JJ S*CW AS 3 *VCMB* t 
2XSTOR1C JJI*CWB5*XSTGR3t JJ) CWB8*XSTGR4i JJ)*CWB9 S+XST0R4OJ?* 


3XST0R5* JJ)*CWB10)*VCMC*12-0*XSTCR2{ J J 3 *(CWC3+XST 0R4( JJ }*CWC4? + 
4XST0RM JJ)*tCWC 5+ XST0R4 ( J J5 *XST0R4{ JJ 1 *CWC6) +XST0R3t J J >* ICWC8* 
5XST0R4t JJ)*CWC9? )) 

PAF 2X2=1 1-O/COEF AS) *t V CMA*12 - 0*XSTOR2 ( J J )#CW A1*XST0R1 t J J )* t 
IC m 34- XST0R4T J J ? *C WA4I 4* XSTOR3 ( JJ) *CWA7 +X ST GR5 { J J 3 * CW A8 I *V CM B* ( 
2XSTGRH JJ3*CWB5+XSTOR3t JJ) * ? CWB8+XSTOR4 t JJ }*CW89 UXSTOR4tJ J 3* 
3XSTOR 5( J J) *C WB 1 0? * VCNC*(2«» □❖XSTGR2 ( JJ ) CWC3*XST 0R4( J J )*CWC4?-i- 
4XST0R It JJ) ❖ ( CWC54- XSTOR4 1 J J) *XST0R4 < JJ ) *CWC& J *XS7QR3 1 J J )* tCWCS* 
5XSTGR4C JJ?*CVJC93 ) ) 

RAF 3X2= t 1-0/CGEFAS3*! VCPA*(2.0*XST0R2 t J J 3 *CW Al+XST OR1 ( J J )* t 
1CWA 3*XSTOR4t JJ? *C WA4? *XSTOR3 t J J3 *CWA7+XSTOR5 t JJ )*CW A8 ) +VCP B* ( 
2XSTOR U J J)*CWB 54- XST0R3 1 J J) *t CWB8 4-XST0R4 { J J 3 * CW B9 ) 4-XST0R4 t J J )* 
3XSTGR5( JJ)*CWB10)*VCPC i M2c0*XSTGR2( J J ) * ( CW C3+XST0R4 ( J J 3* CW C4 ) + 
4XST0R It JJ?*(CWC 5+ XSTGR4C JJ? ❖XSTQR4 CJJ)*CWC6?+XST0R3tJJ J* (CWC8+ 
5XST0R4t JJ)*CWC9? ) ) 

F2=“P( 1)*PAF1X2-P<2> *PAF2X2-Pt33 *PAF3X2 
l-SOMTC l*XST0R2t J J) / 1 X2 WAX**2 ) 

RETURN 

3 PAF iX3= t 1- O/COEF A S3 *1 VCHA^l 2 . O^XSTOR3 (JJI^CWAZ^XSTORliJJ 3* tCWA5 
1+XSTOR 4t J J? WA63+ XSTOR2 ( JJ? #C WA7 +X ST OR5 I J J ) =*CW A9 3+VCMB* (2-0* 

2X S TOR 3( J J I ❖ t C WB 3+ X STGR4 1 J J ) *CWB4) *X$TGR1 1 J J) *1 CW B6+XST0R4UJ 3* 
3XST0R4C J J ? *CWB 73 + XSTGR2 I J J) ❖CCWB8+XSTGR4 «J J) *CW B9)3+VCMC*( 
4XSTQRH JJ?*CliC74-XSTOR2 t JJ) ❖(CWC8*XST0R4 t J J)*CW C93 -t-XST 0R4CJJ 3* 
5XSTOR5S JJ)*CWC103 3 

PAF2X3=iU0/CGEFAS3=M VCNA*t2-0*XSTGR3 (JJ)*CWA2*XSTORHJJ S* ICWA5 
1* XS TOR 4U J 3 *C W A 6 3 * X S TGR2 1 J J ) *C W A7 *X S T BR5 t J J I *CW A9 ) +V CN BM 2 - 0* 
2XSTQR 3t JJ) HZ WB3* XST0R4! JJ) *CWB4 J *XST OR1 i J JJ’M CW B6*XST0R4(J J )❖ 
3XSTGR4C JJ? ❖CWBTJ-t- XSTQR2 { JJ? ❖ ( CWB8-S-XST0R4 ( J J3 *CW B9 ) ) *VCNC* t 
4XST0R It J,H*CWC7*XSTOR2S JJI *(CWC8*XST0R4( JJ3*CWC9?*XSTGR4UJ 3* 
5XSTGR5C JJ3*CWC10) ) 

PAF3X3=t 1-0/CGEFAS? ❖ ( VCPA*t2„0*XSTOR3 t J J S*€W A2*XSTGR 1 1 J J )=MCWA5 
I*XSTGR4i J J 5 *C WA 63 + XSTGR2 i J J) *C WA7 *X$T GR5 ( J J 3 *CW A9 3 *V CP B* < 2 * 0* 

2 XSTDR 3f J J 3 * t C WB 34- X STGR4 1 J JJ ❖CW B4 3 *X ST GR H J J ? * < CW B6*XST GR4 ( J J )* 
3XSTGR4C JJ3 *C WBT) 4- XSTOR2 1 J J? ❖CCWB8+XSTGR4 { J J1 *CW B9 U +VCPC* ( 
4XSTDR 1C JJ? *C WC7*XSTGR2 t J J) * tCWCB*XSTOR4 { J J )*CWC9 HXST0R4U J ?* 


5XST0R 5i JJ3*CWC10) I 

E 2— P < 1 ) *P AF 1 X3-P 12) *P AF2X3-P 1 3 3 *P AF3X 3 
l-$ONTCl*XSTGR3? JJ) /CX3MAX**2) 

RETURN 

4 IH XST0R4C J J) *GE » PEXLWH) GO TO 31 
IF! XST0R4UJT o LE*PEXLMU GO TO 32 

30 PENCO 1=0® 0 
GO TO 33 

31 PENC01=PECXVA 

GO TO I 7501 * 33, 7501 ) „M 

7501 PENCO 1=0® 0 
GO TO 33 

32 PENCO 1=-PECXVA 

GO TO { 7502, 7502, 33) 

7502 PENCO 1=0* 0 

33 DVCMA-2«0 J 5 t XSTGR4{ J J1 #CWAl+4* 0*XST0R4 l JJ) **3*CMA2 

DVCMB=CWB 1+2® G*XST0R4T JJ) *CMB24*3®0*XSTOR4v ? JJ)*XSTGR4UJ )*CMB3 
D VC MC=C HC 1+ 2o 0*X$ T0R4 \ JJS *CMC2 +3 * 0 *XS TO R4 U J ) *XST 0R4 i J J 3* C HC 3 
0 VC NA=C NA 1+ 2* 0 *XSTQR4 l J J) *C NA2 4-3 * 0 *X5 T0R4 < J J I *X ST 0R4 1 J J )£CNA3 
DVCNB=CNB 1+2- 0*XST0R4i JJ) *CNB2 
DVCNC=C NC 1+2* 0*X$TGR4t J J) *CNC2 

D VC P A =C PA 1+ 2® 0* X S T0R4 1 J J) *C P A2 4-3 ® 0 *X S T 0 R4 C J J 1 ** Z* CP A3 
DVCPB=CPB 1+2. G*XSTDR41 J JJ *CPB2 
DVCPC-CPC 1+2* O*XST0R4C JJ3 *CPC2 

VAA-XST0R21 J J] *XSTGR2 f J J? #CWA1 +XSTGR3 C J J 3 *XST 0R3 C J J )*CW A2+XST0R 1 
lUif mST0R2f JjmCWA3+XST0R4f JJ)*CWA4)+XSTGR1{JJ3*XSTGR3UJ) 

2# t C ^WA 5+ X STQR 41 J 3 1 *CViA&)+XST0R2 TJ J) *XST0R3 C JJ)*CWA7+XST0R2{JJ ) 
3*XSTGR5{ JJ1*CWA8+XSTSR3« JJ3*X$TDR5f JJJ*CWA9 
VBB=XSTGRHJJ)*XSTGR1C JJ3*{C«Bl+XST0R4t JJ)*CW82)-i-XST0R3UJ » 
1*XSTGR31 JJ3*ICWB3+XSTOR4C JJ3*CVIB4)+XST0RU JJ5*XSTOR2(JJ)*CWB5 
2+XSTOR If J33 *XSTGR3 C JJS #ICV3B6+XSTGR4 ( JJ) *XSTGR4{ JJ )*CVJB7) +XST0R2 
31 JJ 3*XST0R3C J J) m CWBB+XST0R4 1 3 J 3 *CW B9 3 4-XST0R21 J J )* { XST0R4( JJ )* 
4XST0R5! JJ)*CWB103 + U* JJ) *CWB11 
VCC-XSTORlf JJ3 *XST0R1 { JJ) *( CWCH-XST0R4 ( J J ) *CWC2 ) +X$T0R2( JJ 3 
l^XSTOR 21 JJ) #CC WC3+ XST0R4I JJ) *CWC4 3+XSTGRK JJ)*X$T0R2{JJ »* 


2(CWC5-**XSTQR4UJ)*XST0R4l JJ) *C1C6UXST0R1 ( JJ)*XST0R3( JJ >*CWC7 
34-X5T0R2C J J) *XSTGR3! J J) CWC8+XST0R4 ( J J) *CWC9 H-XST0R3C J J )* { XSTGR4 
4( J J )*XSTGR5f J J) *CWC1 01 *Ut J J) *C WC11 
PAF 1 X4= { lc 0/CDEFAS)*! VAA#DVCMA*VCMA*l XSTORU JJ )*XST0R2<JJ }*CWA4 
1+XSTORli JJ) *XST0R3UJ) *CWA6 ) + VBB*DVCMB-H/ CMB* CXST QR1 U J J*XSTORHJJ S 
2*CWB2*XST0R3l JJ) *XST0R3 I J J) *CW84+XST0R1 < J J ) *XST 0R3 l J J )*2.0 
3*XST0R4UJ)*CWB7+XSTOR2i JJ) *XST0R3 t J J ) *CW B9+X$T0R2{ J J )*XST0R5UJ ) 
4*CWBI0I+VCC*DVCMC + VCMC*< XST0R1 t J J) *XSTORH JJ) *CWC2*XST0R2(JJ ) 
5*XST0R2t JJ) *CWC4+XST0RU JJ) *XSTGR2 ( JJ 1 *2*0*XST0R4( JJ )*CWC6 
6*XST0R2( JJ) #XST0»3 l J J) *CWC9*XST0R3 t J JS *XST0R5 { J J )*CWC103 ? 
7+tVAA*VCMA+VBB*VCWR + VCC*VCMC) *(-l „0 / ( CGEFAS* COEFAS ) J* (4,0* 

8XST0R4! JJ)**?*"Ynu*V0U*CA3*3o0*XST0R4( JJ)*XST0R4( JJ)*CA23+2.0 

9*x$ tor 4< a j s * vm. ; -?c /, ? + c a i ? ) 

PAF 2X4={ l.O/COEFAS) *\ VAA*DVCNA+VCNA*(XSTORU J J)*XST0R2(J J )*CWA4 
UXSTGR l(JJ)*XST0R3iJJ) *CWA6) *VBB*DVCNB+VCNB* (XSTORU JJ)*XSTQR1(JJ ) 
2*CWB2*X$T0R3( JJ)*XST0R3l JJ) *CWB4*XSTGRU J jmST0R3< J J )*2o0 
3*XST0R4C J J) *CV1B7* XST0R2 I J J) SXST0R3 ( J J) *CWB9*XST 0R2 1 J J 5#X$TQR5{ J J ) 
4^CWB 105 -5* VCC VCNC4- VCNC*{ XSTCRI ( J J) #X5T0R1 { JJ )*CW C2+XSTE3R2{ J J ) 
5*XSTGR 2C J J) WC 4* XSTGR1 1 J J) *X STGR2 ( J J ) *2 « 0*X ST GR4 { J J )* Ctf C6 
6*X5TGR2( JJ) *XST0R3( J J) #CWC9+XST0R3 ( J J $ *X$TGR5 { J J )*CWC 10 H 
7^CVAA*VCMA^VBB*VCNB^VCC*VCNC)*(-l o 0/tC0EFAS^CGEFAS ))*(4.0* 

8XSTGR4C J J S **3*YGU*Y0U*CA3*3 o 0*XSTOR4 ( J J 5 *XSTGR4 1 J J )*CA23*2 .0 
9- XS TOR 41 JJ) *YO0*CA2-s-CA12) 

PAF3X4“t 1. 0/C GEFAS)*mA*DVCPA+VCPA*< XSTORU JJ } *XSTGR2( J J )*'CWA4 
l*XSTGRlf JJi *XSTGR3( JJ)*CWA6M-VBB*DVCPB*VCPB*{XST0RUJJ)*XSTORltJJ ) 
2&C W8 2* X 5TGR 3( J J ) * X ST0R3 ( J J 5 *C W B4 *X ST 0 R 1 ( J J ) *X S T OR 3 ( J J ) * 7 . 0 
3&X5TGR 41 JJI^CWB 7 + XSTQR2 ( JJ) *XST0R3 ( JJ I^CV)B9-*-XST 0R2C J J )*XSTGR5(JJ ) 
4*CWB10UVCC*DVCPC*VCPC*( XSTORU JJS*XSTGRUJJ)*CWC2+XST0R2UJ I 
5*XST0R2( JJ1 *C VIC 4*X STORK JJ) *XST0R2* JJ) *2oO*XSTGR4( JJ )*CWC6 
64-XSTOR 21 J J! #*STGR3 ! J J) *CWC9*XST0R3 ( J J)*XST 0R5 ( J J )*CWC1G) ) 
VAA*VCPA-FVBB*VCP3-i-VCG*VCPC) #t^l«.0/ (CGEFAS*COEFAS ) I* (4.0* 

8XSTGR 41 J J ) *#3*YG0* YGU*C A3*3 * 0*X ST0R4 t J J) #XST 0R4 ( J J )* C A23+2 * 0 
9*X$TGR4( JJ) *YGU^C A2-2-CA12 1 
IFCPENCOloGToOo 0) GO TO 9430 
IF( PENCDI* tTo Oo 0) GO TO 9431 


r 



9430 CONTINUE 

EXTPE V-XSTQR4I J J) — PE XLMH 

GO TO 9432 

9431 CONTINUE 

EXTPE V=~ XSTQR4I JJ?*PEXLML 

9432 CONTINUE 

F2--PI l}*PAFlX4-P(2)*PAF2X4-POl*PAF3X4-2*C*PENtGl*EXTPEV-S0MTC0* 
1XST0R4C JJl/tPE XLIM**2? 

RETURN 

5 P AF 1X5= { loO/COEFASS VCMA*I XSTOR2 I JJ) *CW A8*XSTQR3! JJ S*CW A9 H- 
IVCMB^I XSTOR2? JJ? -XSTOR4I JJ? *CWB10H*VCMC# {XSTGR3 UJ >*XST0R4IJ J ) 
2*CWC 101 ? 

PAF 2X5={ UO/COEFAS)*! VCNA*t XSTGR2I JJ) #CU A8*XST0R3I JJ l*CWA9)* 

1 1VCNB*I XSTGR2! JJ) *XSTGR4t JJ? *CWB101+VCNC*IXSTGR3$ JJ ?*XST0R4U J ) 

V 2#CWC 10 I ) 

p PAF3X5=f loO/COEFAS? #S VCPA^f X ST0R2 f JJ? *CW A8*XSTGR3! JJ 3*CWA9K 

XVCPB*C XST0R2C JJ3*X$TOR4I JJ? *CWB10) *VCPC* CXST0R3 1 JJ )#XSTQR4 C J J } 
2*CVJC 101 3 

r E2=-P{13*PAFIX5-P12) *PAF2X5-P<3? *PAF3X5-PI4) 

1- $0MTC2*XST0R5{ JJ? /CX5MAX**2) 

'e RETURN 

END 


h 
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